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Abstract 


Compressive  Spectral  Imaging  is  a  revolutionary  technique  which  senses  the  spatio-spectral  information 
of  a  scene  by  using  2D  coded  projections.  The  underlying  spectral  3D  data  cube  is  then  recovered  using 
compressed  sensing  (CS)  reconstruction  algorithms  which  assume  that  the  underlying  hyperspectral  im¬ 
ages  are  sparse  in  some  representation  basis.  The  great  advantage  of  CSI  is  that  the  required  number  of 
measurements  needed  for  reconstruction  is  far  less  than  that  required  by  traditional  scanning  methods. 
In  practice,  the  Coded  Aperture  Snapshot  Spectral  Imaging  (CASSI)  systems  efficiently  implement  CSI. 
The  compressive  CASSI  measurements  are  often  modeled  as  the  summation  of  coded  and  shifted  versions 
of  the  spectral  voxels  of  the  underlying  scene.  Thus,  each  CASSI  measurement  is  a  highly  structured 
random  projection  of  the  underlying  scene.  The  structure  of  these  projections  is  dictated  by  the  CASSI 
optical  system  whose  only  varying  components  are  the  coded  aperture  entries.  The  coded  apertures 
are  crucial  inasmuch  as  they  determine  the  quality  of  the  image  reconstruction  as  well  as  the  required 
minimum  number  of  FPA  measurements. 

For  almost  a  century,  coded  apertures  have  been  designed  by  the  realization  that  they  are  well  con¬ 
ditioned  when  used  in  least  square  estimation.  In  fact,  commonly  used  coded  apertures  in  CASSI  in¬ 
clude  Hadamard  matrices  Hat  whose  entries  are  (iCv)t?  e  {  — l,l}7VxiV,  Hadamard  S  matrices  Sat  where 
Sat  =  1/2(1  —  Hat),  cyclic  S-matrices  consisting  of  cyclic  permutations  of  a  single  master  codeword,  and 
Bernoulli  random  matrices.  In  this  project,  a  new  generation  of  coded  apertures  for  compressive  spectral 
imaging  are  developed.  The  new  coded  apertures  are  designed  such  that  they  are  not  only  optimal  for 
least  square  methods,  but  also  for  recovery  methods.  We  coin  the  new  class  of  coded  apertures  as  “fl¬ 
ooded  apertures”  since  these  are  optimal  under  the  criteria  often  used  in  the  emerging  techniques  of  CS. 
Notably,  the  t\ -coded  apertures  ensembles  can  be  designed  to  obtain  better  quality  in  the  reconstructed 
images  while  minimizes  the  number  of  required  2D  projections. 

By  designing  the  coded  aperture  patterns  such  that  the  sensing  better  satisfies  the  Restricted  Isome¬ 
try  Property  (RIP)  of  the  CASSI  sensing  matrix.  The  project  thus  establishes  the  CASSI  RIP  constants 
in  terms  of  the  statistical  properties  of  the  coded  aperture  entries  such  that  the  concentration  of  the 
eigenvalues  of  the  sub-matrices  associate  to  the  CASSI  sensing  matrix  is  maximized.  More  specifically, 
the  optimal  i\ -coded  apertures  are  designed  when  their  entries  are  drawn  from:  (a)  Boolean  random 
variables,  (b)  binary  random  variables,  (c)  random  signed  and  unsigned  gray  scale  values,  ^i-coded 
apertures  can  also  be  designed  for  specific  applications  such  as  compressive  spectral  selectivity.  Spec¬ 
tral  imaging  selectivity  is  sought  in  diverse  applications  since  relevant  information  often  lies  within  a 
subset  of  spectral  bands.  Capturing  and  reconstructing  all  the  spectral  bands  in  the  observed  image 
cube,  to  then  throw  away  a  large  portion  of  this  data  is  inefficient.  First,  the  i\ -selective  coded  aperture 
optimization  problem,  in  this  case,  is  shown  to  be  analogous  to  the  optimization  of  a  constrained  mul¬ 
tichannel  filter  bank.  The  optimal  ^i-selective  coded  apertures  allow  the  decomposition  of  the  CASSI 
measurement  into  several  periodic  subsets,  each  having  information  from  only  a  few  selected  spectral 
bands.  This  first  design  approach,  however,  limits  the  selective  spectral  profiles  to  be  periodic  patterns. 
Further,  the  minimum  number  of  shots  is  restricted  to  the  periodicity  of  the  spectral  pattern.  In  most 
practical  applications,  however,  the  spectral  profiles  of  interest  are  not  periodic  and  the  number  of  shots 
is  additionally  restricted  by  the  application  at  hand.  A  second  optimization  approach  for  coded  aperture 
spectral  selectivity  is  then  developed,  where  the  ^i-selective  coded  aperture  design  is  generalized  to  a  more 
general  and  more  effective  mathematical  framework  for  multi-shot  CASSI,  allowing  the  reconstruction  of 
arbitrary  subset  of  bands,  periodic  or  aperiodic,  while  minimizing  the  required  number  of  shots.  The  new 


10 


approach  aims  at  the  optimization  of  ^-selective  coded  aperture  sets  such  that  a  group  of  compressive 
spectral  measurements  is  constructed,  each  with  information  from  a  specific  subset  of  bands.  A  matrix 
representation  of  CAS  SI  is  introduced  permitting  the  optimization  of  spectrally  selective  coded  aperture 
sets.  Further,  each  ^-selective  coded  aperture  set  forms  a  matrix  such  that  rank  minimization  is  used  to 
reduce  the  number  of  CASSI  shots  needed.  Conditions  for  the  code  apertures  are  identified  such  that  a 
Restricted  Isometry  Property  in  the  CASSI  compressive  measurements  is  satisfied  with  higher  probability. 

Finally,  this  project  develops  a  higher  order  precision  model  for  the  optical  sensing  in  CASSI  that 
includes  a  more  accurate  discretization  of  the  underlying  signals,  leading  to  image  reconstructions  less 
dependent  on  calibration.  Further,  the  higher  order  model  results  in  improved  image  quality  reconstruc¬ 
tion  of  the  underlying  scene  than  that  achieved  by  the  traditional  model.  The  proposed  higher  precision 
computational  model  is  also  more  suitable  for  reconfigurable  multi-frame  CASSI  systems  where  multi¬ 
ple  coded  apertures  are  used  sequentially  to  capture  the  hyperspectral  scene.  Several  simulations  and 
experimental  measurements  demonstrate  the  benefits  of  the  new  discretization  model. 
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Introduction 


Compressive  Spectral  Imaging 

Imaging  spectroscopy  involves  the  sensing  of  a  large  amount  of  spatial  information  across  a  multitude 
of  wavelengths.  Traditional  imaging  spectroscopy  sensing  techniques  scan  adjacent  spatial  zones  of  an 
underlying  spectral  scene  and  merge  the  results  to  construct  a  spatio-spectral  data  cube.  Push  broom 
spectral  imaging  sensors,  for  instance,  capture  the  spectral  data  cube  by  using  a  dispersive  element  as  a 
prism  or  grating  and  one  Focal  Plane  Array  (FPA)  measurement.  These  systems  capture  one  snapshot 
per  spatial  line  of  the  scene  [?]  and  the  measures  are  concatenated  to  construct  the  spatio-spectral  data 
cube.  Spectrometers  based  on  optical  band-pass  filters  scan  a  scene  by  tuning  band-pass  filters  in  wave¬ 
length  steps  such  that  a  whole  spectral  region  is  covered  [?]. 

Although  these  traditional  sensing  techniques  can  be  designed  to  measure  the  spectral  information 
with  high  resolution,  they  have  the  disadvantage  that  they  require  to  scan  a  number  of  zones  (spatial 
lines,  spectral  bands)  linearly  in  proportion  to  the  desired  spatial  or  spectral  resolution.  In  contrast, 
Compressive  Spectral  Imaging  (CSI)  senses  2D  coded  projections  of  the  underlying  scene  such  that  the 
number  of  required  projections  is  far  less  than  the  linear  scanning  case.  CSI  exploits  the  fact  that  hyper- 
spectral  images  are  sparse  in  some  basis  representation  [?].  More  formally,  suppose  that  a  hyperspectral 
signal  T  G  M NxNxL ?  Qr  its  vector  representation  /  G  RN'N'L,  is  S  sparse  on  some  basis  such  that 
/  =  ^ 6  can  be  approximated  by  a  linear  combination  of  S  vectors  from  with  S  <C  (N  •  N  •  L).  Here, 
N  x  N  represents  the  spatial  dimensions  and  L  is  the  spectral  depth  of  the  image  data  cube.  CSI  allows  / 
to  be  recovered  from  m  random  projections  with  high  probability  when  m  >  S'log(TV-TV-L)  <C  (N  -N  -L). 

Coded  Aperture  Snapshot  Spectral  Imager  (CAS SI)  is  a  spectral  image  sensor  that  attains  FPA  com¬ 
pressive  measurements  by  using  a  coded  aperture  and  a  dispersive  element.  Indeed,  CAS  SI  has  motivated 
a  great  diversity  of  research  directions  in  areas  such  as  compressive  spectral  classification  [?,  ?],  compres¬ 
sive  tomography  [?],  compressive  holography  [?],  X-ray  compressive  imaging  [?],  compressive  fluorescence 
microscopy,  [?],  spatio-spectral  compressive  super-resolution  [?,  ?],  and  compressive  Raman  spectroscopy 
[?,  ?].  The  projections  measured  in  CASSI  are  given  by  y  =  H /,  where  H  is  an  N(N  +  L  —  1)  x  (TV  -N  -L) 
matrix  whose  structure  is  determined  by  the  coded  apertures  and  the  dispersive  element. 

For  spectrally  rich  scenes  or  very  detailed  spatial  scenes,  a  single  shot  CASSI  measurement  may  not 
provide  a  sufficient  number  of  compressive  measurements.  Further,  increasing  the  number  of  FPA  random 
projections  yields  to  a  less  ill-posed  problem,  however,  each  additional  snapshot  requires  more  integration 
time  at  the  detector.  Each  additional  shot  uses  a  distinct  coded  aperture  that  remains  fixed  during  the 
integration  time  of  the  detector.  This  will  rapidly  increase  the  quality  of  image  reconstruction.  The 
additional  time  required  to  sense  the  supplemental  FPA  projections  can  be  impractical  in  applications 
such  as  hyperspectral  sensing  at  video  rates  [?].  Increasing  the  number  of  FPA  projections  also  reduces 
the  compressive  advantage  of  CSI.  Minimizing  the  required  number  of  FPA  projections  is  thus  a  key  issue 
in  CSI.  The  CASSI  architecture  has  been  recently  modified  to  admit  multiple  snapshots,  each  admitting 
a  different  coded  aperture  pattern.  Multi  shot  CASSI  can  thus  yield  a  less  ill-posed  inverse  problem  and 
consequently  improved  signal  recovering  [?,  ?].  The  time- varying  coded  apertures  can  be  implemented 
using  micro-piezo  motors  [?]  or  through  the  use  of  Digital  Micromirror  Devices  [?,  ?].  Denote  the  ith  FPA 
measurement  as  y£  =  H where  EC  represents  the  effects  of  the  ith  coded  aperture.  The  set  of  K  FPA 
measurements,  each  with  a  different  coded  aperture  is  then  assembled  as  y  =  [(y°)T, . . . ,  (yx_1)T]  • 
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The  projections  in  CASSI  can  be  alternatively  expressed  as  y  =  =  A 0  where  the  matrix  A  =  H\£ 

is  the  sensing  matrix.  The  underlying  data  cube  is  reconstructed  as  /  =  ®  (min#  ||y  —  H\£0||2  +  r||0||i) 
where  H  =  [(H°)T, . . . ,  (HX_1)T]T,  0  is  an  S  sparse  representation  of  /  on  the  basis  and  r  is  a 
regularization  constant. 

Coded  Aperture  Design  in  CASSI 

A  critical  component  in  this  inverse  problem,  is  the  structure  of  the  sensing  matrix  A  =  H\£  as  it  ul¬ 
timately  determines  the  attainable  quality  of  reconstruction.  While  the  optical  architecture  in  CASSI 
partially  imposes  a  well  defined  sparse  structure  to  the  sensing  matrix  A,  the  coded  apertures  used  in 
each  measurement  shot,  ultimately  determine  the  structure  of  A.  The  objective  in  CASSI  is  thus  to 
optimally  design  the  set  of  ^i-coded  apertures  so  as  to  forge  a  structure  on  A  that  minimizes  the  number 
of  FPA  snapshots  while  attaining  the  highest  quality  reconstruction.  To  this  end,  this  project  explores 
the  interplay  between  the  Restricted  Isometry  Property  and  the  set  of  coded  apertures  used  in  CASSI. 

The  Restricted  Isometry  Property  (RIP)  characterizes  the  ” goodness”  of  a  matrix  in  CS,  and  it  is 
used  to  develop  many  theorems  in  CS.  It  establishes  the  conditions  necessary  for  A  such  that  the  t 2 
norm  of  the  underlying  3D  spectral  image  is  approximately  preserved  under  the  transformation  A 0. 
More  precisely,  the  restricted  isometry  constant  Ss  of  the  matrix  A  is  the  smallest  constant  such  that 
(1  —  5s)||0||2  <  ||A0||2  <  1 10| I2 (1  T  [?].  The  RIP  requires  that  all  m  x  \T\  column  submatrices  A|T|  of 
A  are  well  conditioned  for  all  \T\  <  S.  Indeed,  the  RIP  imposes  that  all  the  eigenvalues  of  the  matrices 
Aj^|A|T|  are  in  the  interval  [1  —  Ss  ,  1  +  5S].  The  probability  of  satisfying  this  condition  is  calculated 
by  estimating  the  statistical  distribution  of  the  maximum  eigenvalue  A  max  of  the  matrices  A^  A|T|  —  I 
where  I  is  an  identity  matrix.  The  distribution  of  the  maximum  eigenvalue  Xmax  is  estimated  using 
the  concentration  of  measure  for  random  matrices  developed  in  [?].  The  RIP  condition  also  implies  a 
stable  recovery  of  the  signal  0  from  the  projections  A0  using  an  l\  optimization  algorithm  [?].  Fur¬ 
thermore,  the  RIP  can  be  used  to  determine  bounds  on  the  required  number  of  measurements  needed 
for  successful  CS  reconstruction.  These  bounds  depend  on  the  structure  of  the  underlying  sensing  matrix. 

The  first  approach  to  characterize  the  RIP  in  CASSI  was  given  in  [?],  where  it  was  assumed  that  the 
RIP  for  the  matrix  A  is  satisfied  for  some  constant  Ss ,  then  conditions  on  the  coded  apertures  were  deter¬ 
mined  so  that  the  RIP  is  better  satisfied.  These  results,  however,  do  not  provide  the  explicit  parameters 
for  the  bounds  needed  in  the  RIP,  such  as  the  probability  of  error,  or  the  bound  on  the  minimum  number 
of  FPA  measurements  for  stable  recovery.  Results  in  [?]  also  do  not  exploit  the  RIP  in  the  optimal  design 
of  the  coded  apertures.  In  fact,  commonly  used  coded  apertures  in  CASSI  include  Hadamard  matrices  Hat 
whose  entries  are  ( H^)ij  £  {  —  l,l}iVxiV,  Hadamard  S  matrices  Sat  where  Sat  =  1/2(1  —  Hat)  [?],  cyclic 
S-matrices  consisting  of  cyclic  permutations  of  a  single  master  codeword  [?],  and  Bernoulli  random  ma¬ 
trices  [?,  ?].  The  use  of  these  coded  apertures  in  CASSI  has  been  principally  motivated  by  the  realization 
that  they  are  well  conditioned  when  used  in  least  square  estimation  [?,  ?].  However,  these  code  designs 
do  not  fully  exploit  the  rich  theory  of  CS.  In  particular,  they  do  not  exploit  the  RIP  condition  or  the 
concentration  of  measure  of  the  respective  random  submatrices  of  A  to  define  optimal  coded  aperture  sets. 

The  set  of  optimal  coded  apertures  under  the  criterion  of  euclidian  and  i\  distances  are  coined  here  as 
£1  coded  apertures.  There  are  several  types  of  t\  codes,  boolean  (  Hij  E  {0, 1}),  binary  (. Hij  E  {  —  1, 1}), 
grayscale  G  { 7^4,  •  •  • ,  1}  and  G  {-1,  ,  •  •  • ,  1}),  and  other  possible 

combinations.  The  optimal  ^i-coded  apertures  are  also  designed  for  spectral  selectivity.  Spectral  se¬ 
lectivity  is  of  interest  in  many  applications,  including  wide-area  airborne  surveillance,  remote  sensing, 
and  tissue  spectroscopy  in  medicine.  The  optimal  spectral  bands  in  airborne  surveillance,  for  instance, 
depend  on  atmospheric  conditions,  time  of  day,  the  targets  of  interest,  and  the  background  against  which 
the  targets  are  viewed  [?].  In  these  applications,  the  spectral  signatures  of  interest  live  in  a  spectral  band 
subspace.  Efforts  placed  on  acquiring  the  entire  spectral  image  cube,  to  then  throw  away  a  large  portion 
of  this  data  is  wasteful  in  many  regards. 

A  first  approach  for  spectral  selectivity  is  developed  when  the  selective  spectral  profiles  are  limited 
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to  periodic  patterns  and  the  minimum  number  of  shots  is  restricted  to  the  periodicity  of  the  spectral 
pattern.  In  most  practical  applications,  however,  the  spectral  profiles  of  interest  are  not  periodic  and  the 
number  of  shots  is  restricted  by  the  application.  Then  a  second  approach  is  developed  which  establishes 
a  more  general  and  more  effective  mathematical  framework  for  multi-shot  CAS  SI  and  the  corresponding 
algorithms  for  i\ -selective  coded  aperture  optimization  that  allow  the  reconstruction  of  arbitrary  subset 
of  bands,  periodic  or  aperiodic,  whilst  minimizing  the  required  number  of  shots. 

The  organization  of  the  document  is  as  follows: 

Chapter  1  describes  the  principal  theoretical  concepts  involved  in  the  project.  The  concept  of  spar¬ 
sity  is  illustrated  and  the  physical  phenomenon  in  CASSI  is  presented.  The  techniques  to  recover  the 
underlying  signal  and  the  Restricted  Isometry  Property  are  also  presented. 

Chapter  2  presents  the  coded  aperture  design  for  spectral  selectivity  in  CASSI.  In  many  applications 
selective  spectral  imaging  is  sought  since  relevant  information  often  lies  within  a  subset  of  spectral  bands. 
Capturing  and  reconstructing  all  the  spectral  bands  in  the  observed  image  cube,  to  then  throw  away 
a  large  portion  of  this  data  is  inefficient.  To  this  end,  this  chapter  extends  the  concept  of  CASSI  to  a 
system  admitting  multiple  shot  measurements  which  leads  not  only  to  higher  quality  of  reconstruction, 
but  also  to  spectrally  selective  imaging  when  the  sequence  of  code  aperture  patterns  is  optimized.  The 
coded  aperture  optimization  problem  is  shown  to  be  analogous  to  the  optimization  of  a  constrained  mul¬ 
tichannel  filter  bank.  The  optimal  coded  apertures  allow  the  decomposition  of  the  CASSI  measurement 
into  several  subsets,  each  having  information  from  only  a  few  periodical  selected  spectral  bands.  The  rich 
theory  of  compressive  sensing  is  used  to  effectively  reconstruct  the  spectral  bands  of  interest  from  the 
measurements.  A  number  of  simulations  are  developed  to  illustrate  the  spectral  imaging  characteristics 
attained  by  optimal  coded  aperture.  This  chapter  limits  the  spectral  desired  bands  to  periodical  patterns, 
one  more  generalized  method  is  presented  in  Chapter  3. 

Chapter  3  extends  the  selective  coded  apertures  developed  in  Chapter  2  to  admit  periodical  and 
non-periodical  desired  spectral  profiles.  More  specifically,  a  new  coded  aperture  design  framework  for 
multi-frame  Code  Aperture  Snapshot  Spectral  Imaging  (CASSI)  system  is  presented.  It  aims  at  the  op¬ 
timization  of  coded  aperture  sets  such  that  a  group  of  compressive  spectral  measurements  is  constructed, 
each  with  information  from  a  specific  subset  of  bands.  A  matrix  representation  of  CASSI  is  introduced 
permitting  the  optimization  of  spectrally  selective  coded  aperture  sets.  Further,  each  coded  aperture 
set  forms  a  matrix  such  that  rank  minimization  is  used  to  reduce  the  number  of  CASSI  shots  needed. 
Conditions  for  the  coded  apertures  are  identified  such  that  a  Restricted  Isometry  Property  in  the  CASSI 
compressive  measurements  is  satisfied  with  higher  probability.  Simulations  show  higher  quality  of  spec¬ 
tral  image  reconstruction  than  those  attained  by  systems  using  Hadamard  or  random  coded  aperture  sets. 

Chapter  6  presents  a  more  precise  model  for  the  phenomenon  in  CASSI.  The  compressive  CASSI 
measurements  are  often  modeled  as  the  summation  of  coded  and  shifted  versions  of  the  spectral  voxels 
of  the  underlying  scene.  This  coarse  approximation  of  the  analog  CASSI  sensing  phenomena  is  then 
compensated  by  calibration  preprocessing  prior  to  signal  reconstruction.  This  Chapter  develops  a  high¬ 
er  order  precision  model  for  the  optical  sensing  phenomenon  in  CASSI  that  includes  a  more  accurate 
discretization  of  the  underlying  signals,  leading  to  image  reconstructions  less  dependent  on  calibration. 
Further,  the  higher  order  model  results  in  improved  image  quality  reconstruction  of  the  underlying  scene 
than  that  achieved  by  the  traditional  model.  The  proposed  higher  precision  computational  model  is  also 
more  suitable  for  reconfigurable  multi-frame  CASSI  systems  where  multiple  coded  apertures  are  used 
sequentially  to  capture  the  hyperspectral  scene.  Several  simulations  and  experimental  measurements 
demonstrate  the  benefits  of  the  new  discretization  model. 

Chapter  5  discusses  the  implementation  of  the  CS  theory  in  building  a  testbed  optical  imaging  system. 
This  testbed  is  suitable  for  hyperspectral  imaging  applications  wherein  both  the  spatial  and  spectral 
information  of  the  imaging  scene  are  of  interest.  We  demonstrate  the  feasibility  of  this  testbed  by 
developing  a  Digital-Micromirror-Device-based  Snapshot  Spectral  Imaging  (DMD-SSI)  system,  which 
implements  CS  measurement  processes  on  the  3D  spatial/spectral  information  of  the  imaging  scene. 
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Tens  of  spectral  images  can  be  reconstructed  from  the  DMD-SSI  system  simultaneously  without  any 
mechanical  or  temporal  scanning  processes. 


Impact  of  Research  and  Contribution 

The  contributions  of  the  project  are  mainly  related  with  the  optimal  design  of  coded  apertures  in  com¬ 
pressive  spectral  imaging  and  its  corresponding  testbed  implementation.  The  innovation  of  the  code 
design  presented  here  is  the  design  of  the  coded  aperture  based  not  only  the  1 2  norm  but  also  based  on 
the  £1  distance.  The  coded  apertures  ensembles  are  calculated  to  obtain  spectral  selectivity,  minimum 
RIP  constants,  minimum  number  of  shots.  Further  the  codes  are  also  designed  to  be  used  in  applications 
as  compressive  classification,  spatial  and  spectral  super-resolution.  The  optimization  of  the  coded  aper¬ 
tures  consists  on  to  determine  the  spatial  structure  of  the  patterns  to  obtain  the  maximum  quality  in  the 
reconstructed  images  while  it  uses  the  minimum  number  of  measurements.  Additionally,  a  more  precise 
discretization  model  for  CASSI  is  developed.  The  new  model  is  more  suitable  for  the  calibration  of  the 
coded  apertures  in  multi-frame  CASSI. 
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Chapter  1 


Compressive  Spectral  Imaging 


1.1  Sparsity 

Compressive  sensing  exploits  the  fact  that  most  of  the  signals  are  sparse  in  some  representation  basis 
Hyperspectral  images  are  sparsified  using  a  separate  basis  for  the  spatial  axes  ^2 d  and  other  for 
the  spectral  axis  ^id  [?]•  Commonly,  the  2D  wavelet  transform  is  used  to  sparsify  the  images  at  a 
specific  wavelength.  The  spectral  axis  can  be  sparsified  using  a  ID  Wavelength  Transform  or  the  Cosine 
Transform.  Thus,  the  hyperspectral  image  f  is  expressed  as  f  =  ^2 d  ®  dO  where  0  is  the  Kronecker 
product.  Other  basis  can  also  be  used  as  the  Principal  Component  Vectors  or  pre-trained  dictionaries, 
however,  these  representations  require  previous  knowledge  about  the  underlying  signals.  Figures  1.1 
and  1.2  illustrate  the  sparse  representation  coefficients  for  three  different  basis  representation  and  two 
different  data  bases.  The  first  representation  (Figures  1.1(b)  and  1.2(b))  uses  a  ID  Wavelet  Transform  to 
represent  the  data  cubes,  thus  f  =  ^idO  .  The  second  case  (Figures  1.1(c)  and  1.2(c))  uses  a  2D  Wavelet 
Transform  is  used  such  that  f  =  ^2 dO-  The  third  approach  uses  the  Kronecker  of  the  2D  wavelet  basis 
and  the  Cosine  Transform.  It  can  be  observed  clearly  that  the  Kronecker  basis  produces  the  most  sparse 
representation  (Figures  1.1(d)  and  1.2(d)). 


(a)  Original  datacube 


(b)  ID  Wavelet  basis 


(c)  2D  Wavelet  basis 


(d)  Kronecker  basis 


Figure  1.1:  Compression  basis  comparison,  (b)  ID  Wavelet,  (c)  2D  Wavelet,  and  (d)  Kronecker  (DCT-2D 
Wavelet)  basis  representation. 


To  illustrate  the  effect  of  the  sparse  basis  representation,  Figures  1.4  -  1.10  depict  approximations  of 
hyperspectral  images  using  different  levels  of  sparsity.  These  approximations  are  calculated  by  expressing 
the  underlying  signal  in  a  determined  representation  basis,  then  the  estimated  coefficients  are  sorted  by 
magnitude.  A  given  percentage  of  the  sorted  coefficients  are  held  and  the  remanning  are  set  to  zero. 
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a)  Original  datacube 

(b)  ID  Wavelet  basis 

(c)  2D  Wavelet  basis 

(d)  Kronecker  basis 

Figure  1.2:  Compression  basis  comparison,  (a)  Original  scene,  (b)  ID  Wavelet,  (c)  2D  Wavelet,  (d) 
Kronecker  (DCT-2D  Wavelet) 


The  resulting  coefficients  are  transformed  back  to  the  original  space  to  generate  the  approximate  signal. 
Specifically,  Figures  1.3  and  1.7  depicts  two  data  cubes  used  to  illustrate  the  sparsity  concepts.  Figures 
1.4-Figures  1.6  illustrate  the  results  for  the  first  data  cube,  and  Figures  1.8-Figures  1.10  shows  the  results 
for  the  second  data  cube.  The  sparsity  results  are  shown  for  the  spectral  bands  1,  3,  5  and  7  of  the 
original  data  cubes. 

Figures  1.4(b)-(d)  illustrate  the  approximations  by  holding  only  the  1,  5,  and  10  percent  of  the  samples 
of  the  underlying  signals  in  the  first  data  cube.  Figures  1.4(b)-(d)  use  the  ID  Wavelet  as  representation 
function.  Figures  1.5(b)-(d)  and  Figures  1.6(b)-(d)  illustrate  the  results  for  the  2D  Wavelet  and  Kronecker 
representation,  respectively.  Figures  1.8(b)-(d)  illustrate  the  approximations  by  holding  only  the  1,  5, 
and  10  percent  of  the  samples  of  the  signals  in  the  second  data  cube.  Again,  Figures  1.8(b)-(d)  use  the  ID 
Wavelet  as  representation  function.  Figures  1.9(b)-(d)  and  Figures  1.10(b)-(d)  illustrate  the  results  for 
the  2D  Wavelet  and  Kronecker  representation,  respectively.  Notice  as  the  images  can  be  approximated 
with  high  quality  using  only  a  small  percentage  of  the  original  number  of  samples.  Compressive  spectral 
imaging  takes  advantages  of  this  sparsity  phenomenon  to  solve  the  inverse  problem  of  recovering  a  given 
signal  from  its  random  projections. 

1.2  Sensing  Problem 

The  sensing  physical  phenomena  in  CASSI  is  strikingly  simple,  yet  it  adheres  to  the  incoherence  principles 
required  in  CS.  In  its  simplest  form,  CASSI  measurements  are  realized  optically  by  a  coded  aperture,  a 
dispersive  element  such  as  a  prism,  and  a  focal  plane  array  detector  [?,  ?].  The  coding  is  applied  to  the 
(spatial-spectral)  image  source  density  fo(x;  y]  A)  by  means  of  a  coded  aperture  T(x]y)  as  realized  by  the 
CASSI  system  depicted  in  Fig.  1.11,  where  (x;  y )  are  the  spatial  coordinates  and  A  is  the  wavelength  [?]. 
The  resulting  coded  field  fi(x;y;X)  is  subsequently  modified  by  a  dispersive  element  before  it  impinges 
onto  the  FPA  detector.  The  compressive  measurements  across  the  FPA  are  realized  by  the  integration 
of  the  dispersed  field  f2(x;y;  A)  over  the  detector’s  spectral  range  sensitivity. 

It  was  shown  in  [?,  ?]  that  the  discretized  output  at  the  detector  follows 

Yji  =  2/c=0  J~j(£+k)(k)Tj(£+k)  +  Mji  (1-1) 

where  Yj#  is  the  intensity  measured  at  the  j,  i  position  of  the  detector  whose  dimensions  are  Nx  (N+L— 1), 
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Figure  1.3:  Original  hyperspectral  signal.  The  first  row  contains  the  spectral  bands  1-  4,  and  the  second 
row  contains  the  spectral  bands  5-8. 


T  is  the  discretized  spectral  data  cube,  L  is  the  number  of  spectral  bands,  Tji  is  the  coded  aperture, 
and  l is  the  noise  of  the  system.  Assume  that  the  band-pass  filter  of  the  instrument  limits  the  spectral 
components  between  Ai  and  A2.  If  the  pixel  width  of  the  detector  and  of  the  coded  aperture  are  both 
equal  to  A,  then  the  number  of  resolvable  bands  L  is  limited  by  L  —  aX2^Xl  where  aX  is  the  dispersion 
induced  by  the  prism.  The  spectral  resolution  is  limited  by  The  horizontal  and  vertical  spatial  reso¬ 
lutions  are  limited  by  A  and  the  number  of  spatially  resolvable  pixels  of  the  underlying  scene  is  N  x  N. 

The  sensing  mechanism  is  illustrated  in  Fig.  1.12  and  further  described  as  follows.  Since  the  spectral 
dispersion  in  the  prism  affects  the  field  along  one  spatial  dimension,  a  slice  CAS  SI  model  is  sufficient  to 
characterize  the  sensor.  Figure  1.12  shows  a  slice  F  of  the  data  cube  T  along  the  spectral  plane  and 
the  qth  spatial  dimension  of  dispersion.  In  this  case  (1.1)  would  represent  the  model  of  a  slice  of  CAS  SI. 
The  slice  is  first  coded  in  amplitude  by  the  coded  aperture  T  where  the  code  is  “block” or  “unblock”.  In 
essence,  when  a  “block”  code  aperture  element  is  encountered,  the  energy  along  the  entire  row  of  the  slice 
is  “punched  out” .  The  coded  field  is  then  sheared  along  the  spectral  dimension  as  it  transverse  the  prism 
as  illustrated  in  Fig.  1.12.  Finally,  the  coded  and  dispersed  field  is  “collapsed”  in  the  spectral  dimension 
by  integration  in  the  FPA  detector.  Notice  how  the  various  optical  phenomena  realize  compressive  linear 
projections  onto  the  detector.  To  attain  the  2D  CASSI  sensor  model,  one  must  then  replicate  the  above 
slice  model  for  all  the  remaining  slices  of  the  cube,  to  a  complete  2D  detector. 

For  spectrally  rich  scenes  or  very  detailed  spatial  scenes,  a  single  shot  FPA  measurement  is  not 
sufficient  and  additional  shots  are  required,  each  with  a  distinct  code  aperture  that  remains  fixed  during 
the  integration  time  of  the  detector.  Time-varying  coded  apertures  can  be  realized  by  a  spatial  light 
modulator  (SLM)  or  by  a  lithographic  mask  actioned  by  a  piezoelectric  device  [?,?,?].  It  was  also  shown 
in  [?,  ?]  that  the  ensemble  of  say  K  «  L  FPA  shots  in  1-D  vectorized  form  y  =  [yjf, . . .  ,y^_J  can 
be  rewritten  in  the  standard  form  of  an  underdetermined  system  of  linear  equations 

y  =  A0  =  H\I>0  +  LJ  (1.2) 

where  A  =  IFF  e  M^xLiv2  is  the  CASSI  sensing  matrix,  0  is  a  sparse  representation  of  the  data 
cube  in  a  3-dimensional  base  \£,  and  uj  represents  the  noise  of  the  system.  The  matrix  H  in  (1.2) 
accounts  for  the  effects  of  the  coded  aperture  and  the  prism.  The  sensing  matrix  A  thus  couples  H 
with  the  representation  basis  9.  The  coded  aperture  is  considered  binary  and  the  dispersive  element  is 
considered  linear.  In  practice,  it  is  necessary  take  into  account  the  various  optical  artifacts  and  non-ideal 
characteristic  of  the  optical  system.  The  underlying  computational  concepts  are  general  and  are  thus,  in 
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(c)  W1D  (5%)  (d)  W1D  (10%) 


Figure  1.4:  Sparse  hyperspectral  signal  representation  using  the  ID  Wavelet  basis  function.  The  spectral 
bands  in  (a)  are  represented  using  (b)  1%  (c)  5%,  and  (d)  10%  of  the  sparse  representation  coefficients. 
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(a)  Original  spectral  bands 


(b)  W2D  (1%) 


(c)  W2D  (5%)  (d)  W2D  (10%) 


Figure  1.5:  Sparse  hyperspectral  signal  representation  using  the  2D  Wavelet  basis  function.  The  spectral 
bands  in  (a)  are  represented  using  (b)  1%  (c)  5%,  and  (d)  10%  of  the  sparse  representation  coefficients. 


(a)  Original  spectral  bands 


(b)  Kronecker  (1%) 


(c)  Kronecker  (5%)  (d)  Kronecker  (10%) 


Figure  1.6:  Sparse  hyperspectral  signal  representation  using  the  Kronecker  basis  function  (DCT-2D 
Wavelet).  The  spectral  bands  in  (a)  are  represented  using  (b)  1%  (c)  5%,  and  (d)  10%  of  the  sparse 
representation  coefficients. 


(a)  Original  datacube 


Figure  1.7:  Original  hyperspectral  signal.  The  first  row  contains  the  spectral  bands  1-  4,  and  the  second 
row  contains  the  spectral  bands  5-8. 


principle,  applicable  to  imaging  with  FPAs  sensitive  to  any  region  of  the  visible  and  IR  electromagnetic 
spectrum. 


1 . 3  Reconstruction 

The  underlying  hyperspectral  signal  f  G  'Rn‘2l  must  be  recovered  from  the  set  of  random  projections 
y  =  Af  G  Mm  where  m  «  N2L.  Thus,  the  signal  recovery  in  CAS  SI  entails  solving  an  underdetermined 
linear  system  of  equations.  More  specifically,  the  signal  f  is  estimated  as 


f  =  4/  (  min  |  |y  —  H^l 

6 


+  t\\0\ 


(1.3) 


where  H  =  [Hq\  . . . ,  H^_1]  ,  6  is  an  S  sparse  representation  of  f  on  the  basis  and  r  is  a  regularization 
constant.  Solving  (1.3)  relies  on  the  rich  theory  of  CS  which  offers  a  number  of  efficient  reconstruction  al¬ 
gorithms,  including  Matching  Pursuit  type  algorithms  [?]  or  interior-point  methods  such  as  GPSR  [?,  ?]. 
In  the  GPSR  algorithm,  for  example,  each  iteration  requires  in  the  order  of  0(KNAL)  operations.  Hence, 
the  computational  complexity  scales  rapidly  with  N  and  L.  The  probability  of  correct  reconstruction 
of  the  signal  f  is  completely  determined  by  the  structure  of  the  sensing  matrix  A  =  H\£  G  xLN 

where  N 2  is  the  size  of  the  FPA  sensor,  L  is  the  number  of  spectral  bands,  and  K  is  the  number  of 
measurement  shots,  with  K  «  L.  Increasing  the  number  of  FPA  random  projections  yields  to  a  less 
ill-posed  problem,  however,  each  additional  snapshot  requires  more  integration  time  at  the  detector.  The 
additional  time  required  to  sense  the  supplemental  FPA  projections  can  be  prohibited  in  applications 
such  as  hyperspectral  sensing  at  video  rates  [?].  Increasing  the  number  of  FPA  projections  also  reduces 
the  compression  advantage  of  CSI.  Minimizing  the  required  number  of  FPA  projections  is  thus  of  inter¬ 
est  in  CSI  applications.  The  required  number  of  projections  in  CASSI  to  correct  image  reconstruction 
is  determined  by  the  Restricted  Isometry  Property  (RIP)  of  the  matrix  A.  The  RIP  property  can  be 
modified  by  designing  the  set  of  coded  aperture  patterns.  The  optimal  set  of  code  aperture  patterns  is 
designed  such  that  the  RIP  property  is  satisfied  using  the  minimum  number  of  patterns. 


A  more  appropriate  alternative  to  increase  the  probability  of  correct  reconstruction  is  by  designing 


22 


(a)  Original  spectral  bands 


(b)  W1D  (1%) 


(c)  W1D  (5%) 


(d)  W1D  (10%) 


Figure  1.8:  Sparse  hyperspectral  signal  representation  using  the  ID  Wavelet  function.  The  spectral  bands 
in  (a)  are  represented  using  (b)  1%  (c)  5%,  and  (d)  10%  of  the  sparse  representation  coefficients. 
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(c)  W2D  (5%)  (d)  W2D  (10%) 


Figure  1.9:  Sparse  hyperspectral  signal  representation  using  the  2D  Wavelet  function.  The  spectral  bands 
in  (a)  are  represented  using  (b)  1%  (c)  5%,  and  (d)  10%  of  the  sparse  representation  coefficients. 
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(c)  Kronecker  (5%)  (d)  Kronecker  (10%) 


Figure  1.10:  Sparse  hyperspectral  signal  representation  using  the  Kronecker  (DCT-2D  Wavelet)  function. 
The  spectral  bands  in  (a)  are  represented  using  (b)  1%  (c)  5%,  and  (d)  10%  of  the  sparse  representation 
coefficients. 
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Figure  1.11:  Compressive  CAS  SI  sensor  components  developed  at  Duke  University. 


Figure  1.12:  Illustration  of  the  spectral  data  flow  in  CASSI.  The  qth  slice  of  the  data  cube  T  with  11 
spectral  components  is  coded  by  a  row  of  the  code  aperture  t  and  dispersed  by  the  prism.  The  detector 
captures  the  intensity  y  by  integrating  the  coded  light. 


the  structure  of  the  sensing  matrix  A  =  H\I/  as  it  ultimately  determines  the  attainable  quality  of  re¬ 
construction.  While  the  optical  architecture  in  CASSI  partially  imposes  a  well  defined  sparse  structure 
to  the  sensing  matrix  A,  the  coded  apertures  used  in  each  measurement  shot,  ultimately  determine  the 
structure  of  A.  The  objective  in  CASSI  is  thus  to  optimally  design  the  set  of  coded  apertures  so  as  to 
forge  a  structure  on  A  that  minimizes  the  number  of  FPA  snapshots  while  attaining  the  highest  quality 
reconstruction.  To  this  end,  this  project  explores  the  interplay  between  the  Restricted  Isometry  Property 
and  the  set  of  coded  apertures  used  in  CASSI. 
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Chapter  2 


Coded  Aperture  Optimization  for 
Spectrally  Agile  Compressive 
Imaging 

2.1  Introduction 

Coded  Aperture  Snapshot  Spectral  Imaging  (CASSI),  first  introduced  in  [?],  is  a  remarkable  imaging 
architecture  that  allows  capturing  spectral  imaging  information  of  a  3D  cube  with  just  a  single  2D  mea¬ 
surement  of  the  coded  and  spectrally  dispersed  source  field.  The  CASSI  architecture  has  been  studied 
extensively  in  [?,  ?].  It  turns  out  that  the  coded  measurements  are  mathematically  equivalent  to  com¬ 
pressive  projections  which  are  at  the  heart  of  the  emerging  field  of  compressive  sensing  (CS).  In  CS, 
traditional  sampling  is  replaced  by  measurements  of  inner  products  with  random  vectors.  In  CASSI,  the 
random  projections  are  equivalent  to  spectrally  dispersed  coded  fields  that  are  integrated  (projected)  by 
array  detectors  with  broad  spectral  response.  The  spectral  images  are  then  reconstructed  by  solving  an 
inverse  problem  such  as  a  linear  program  or  a  greedy  pursuit  in  a  basis  where  the  under  sampled  signals 
admit  sparse  representations.  The  key  idea  in  CS  reconstruction  is  the  realization  that  most  signals 
encountered  in  practice  are  sparse  in  some  domain  and  the  theory  of  CS  exploits  such  sparsity  to  dictate 
that  far  fewer  sampling  resources  than  traditional  approaches  are  needed  [?,  ?,  ?]. 

Spectral  image  cubes  are  particularly  well  suited  for  sparse  representation  as  images  across  different 
wavelengths  exhibit  strong  correlation  [?,  ?,  ?].  More  formally,  suppose  a  hyperspectral  signal  F  G 
'j^NxMxL ^  Qr  -^s  vector  representation  f  G  7 ZN'M’L,  is  S  sparse  on  some  basis  T  =  Ti(S)T2(S)T3,  such 
that  f  =  T0  can  be  approximated  by  a  linear  combination  of  S  vectors  from  T  with  S  <C  (N  •  M  •  L). 
The  operator  ®  is  the  Kronecker  product  and  T  is  the  Kronecker  basis  representation  of  f  [?] .  N  x  M 
represents  the  spatial  dimensions  and  L  is  the  spectral  depth  of  the  image  cube.  The  theory  of  compressive 
sensing  shows  that  f  can  be  recovered  from  m  random  projections  with  high  probability  when  m  < 
5  log  (AT  •  M  •  L)  <  (TV  •  M  •  L).  The  dimension  of  the  Kronecker  basis  set  T  is  7 zN'M'LxN'M'L  sucj1 
that  Ti  G  7 ZNxN ,  T2  G  TZMxM  and  T3  G  lZLxL .  The  projections  are  given  by  g  =  H where 
H  is  an  N(M  +  L  —  1)  x  (N  •  M  •  L)  random  measurement  matrix  with  its  rows  incoherent  with  the 
columns  of  T.  Commonly  used  random  measurement  matrices  for  CS  are  random  Gaussian  matrices 
(H ij  G  {A/"(0, 1/n)}),  Rademacher  matrices  (H^-  G  {Tl/n1/2})  and  partial  Fourier  matrices.  The  random 
projection  matrices  in  CASSI  are  realized  by  the  coded  aperture  and  the  dispersive  element.  Figure  2.1 
illustrates  the  principles  behind  CASSI,  where  a  6  x  6  x  8  spectral  data  cube  with  16  non  zero  spectral 
components  is  depicted.  The  spectral  components  are  coded  by  a  coded  aperture  where  the  white  pixels 
block  the  impinging  light  and  black  pixels  permit  light  to  pass  through.  A  prism  disperses  the  coded 
light  and  a  detector  integrates  the  intensity  of  the  wave  light.  The  pixel  measurements  at  the  detector 
are  proportional  to  the  linear  combinations  of  the  coded  spectral  components. 

This  chapter  extends  the  CASSI  spectral  imaging  concept  to  a  spectral  imaging  architecture  admitting 
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Figure  2.1:  The  principles  behind  CASSI  imaging.  A  6  x  6  x  8  spectral  data  cube  with  16  non  zero 
spectral  components  is  coded  by  the  coded  aperture  and  dispersed  by  the  prism.  The  detector  integrates 
the  intensity  of  the  resulting  light  wave.  Each  pixel  at  the  detector  contains  a  coded  linear  combination 
of  the  spectral  information  from  the  respective  data  cube  slice. 


multiple  measurement  shots.  The  multiple  measurements  are  attained  as  separate  FPA  measurements, 
each  with  a  distinct  coded  aperture  that  remains  fixed  during  the  integration  time  of  the  detector.  There 
are  several  advantages  to  multiple  shots  [?].  First,  the  number  of  compressive  measurements  in  CASSI 
may  not  meet  the  minimum  needed  for  adequate  reconstruction.  Compressive  sensing  dictates  that  the 
number  of  measurements  must  be  in  excess  of  S' log  (A"  •  M  •  L).  Failure  to  collect  a  sufficient  number  of 
measurements  leads  to  a  severely  underdetermined  inverse  problem  and  to  inadequate  signal  reconstruc¬ 
tion.  With  each  shot,  CASSI  invariably  collects  N  •  M  additional  measurements  regardless  of  the  spectral 
depth  and  sparsity  of  the  source.  For  spectrally  rich  scenes  or  very  detailed  spatial  scenes,  a  single  shot 
CASSI  measurement  may  not  provide  a  sufficient  number  of  compressive  measurements.  Increasing  the 
number  of  measurement  shots  will  multiply  the  number  of  measurements,  thus  rapidly  overcoming  such 
limitations. 

A  second  advantage  to  multiple  shots  is  that  spectral  selectivity  can  be  attained  by  coded  aperture 
design.  Notably,  the  coded  aperture  patterns  can  be  designed  so  as  to  maximize  the  information  content 
on  a  pre-specified  subset  of  spectral  bands  of  particular  interest.  Spectral  selectivity  is  a  characteristic 
of  interest  in  many  applications,  including  wide-area  airborne  surveillance,  remote  sensing,  and  tissue 
spectroscopy  in  medicine.  The  optimal  spectral  bands  in  airborne  surveillance,  for  instance,  depend  on 
atmospheric  conditions,  time  of  day,  the  targets  of  interest,  and  the  background  against  which  the  targets 
are  viewed.  The  development  of  dynamically  programable  multi-spectral  imaging  sensors  are  of  signifi¬ 
cant  interest,  particularly  for  their  use  in  small  unmanned  aerial  vehicles.  Similarly,  spectrally  selective 
imaging  is  becoming  widely  used  in  medicine.  Spectral  imaging  of  the  retina,  for  instance,  offers  a  route  to 
non-invasive  characterization  of  the  biochemical  and  metabolic  processes  within  the  retina  retinopathies 
[?].  In  the  applications  described  above,  and  in  many  other  applications,  the  spectral  signatures  of  interest 
live  in  a  spectral  band  subspace.  Efforts  placed  on  acquiring  the  entire  spectral  image  cube,  to  then  throw 
away  a  large  portion  of  this  data  is  wasteful  in  many  regards.  This  Chapter  focuses  on  overcoming  these 
inefficiencies,  by  developing  an  optimization  framework  for  the  design  of  coded  apertures  for  spectrally 
selective  imaging. 

In  order  to  fully  exploit  the  advantages  of  multishot  coded  apertures,  this  Chapter  also  develops 
a  new  model  describing  the  coded  aperture  optical  system.  The  original  CASSI  model  describes  the 
output  of  the  system  as  a  matrix  operation  acting  on  the  input  source  field,  where  all  the  optical  effects 
are  implicitly  embedded  in  the  matrix  operation.  The  CASSI  model  derived  here  describes  the  output 
of  the  system  as  a  function  of  a  “code-aperture  vector”  that  explicitly  acts  on  a  “reordered”  input 
source.  In  particular,  a  vector  representation  of  a  slice  of  the  data  cube  is  re-ordered  into  a  matrix 
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Figure  2.2:  The  CAASI  system  seen  as  a  multi-channel  filter  bank.  A  set  of  measurements  {gz}  are 
captured  with  a  corresponding  set  of  optimal  aperture  codes  {w1, . . . ,  wK}.  The  mapping  B  reorders 
the  K  sets  of  measurements  {g2}  into  the  V  sets  { g 2}.  Each  subset  {g1}  is  sufficient  to  independently 
reconstruct  the  desired  spectral  band  subset  {FaJ. 

structure,  implicitly  accounting  for  the  spectral  shearing  operator  of  the  CASSI  system.  The  reordering 
operator,  in  turn,  decouples  the  code-aperture  effects  of  CASSI  from  the  spectral  shearing,  thus  allowing 
the  representation  of  the  output  of  the  CASSI  system  as  a  simple  matrix  multiplication  of  the  reordered 
input  data  with  a  vector  whose  elements  are  the  coded  aperture  variables  to  be  optimized.  In  this  manner, 
the  coded  aperture  optimization  problem  is  shown  to  be  analogous  to  a  multi-channel  digital  filter  bank 
optimization  problem  where  the  filter  coefficients  represent  the  aperture  codes  and,  consequently,  the 
optimization  is  constrained  to  produce  binary- valued  outputs.  The  solution  to  the  optimization  problem 
at  hand  is  made  feasible  by  the  new  CASSI  system  representation.  Using  the  resultant  optimal  coefficients, 
a  filter  bank  decomposition  of  the  data  cube  can  be  realized  such  that  it  is  possible  to  assemble  a  reduced 
set  of  compressive  measurements  which  are  sufficient  to  reconstruct  a  desired  subset  of  spectral  bands, 
and  at  the  same  time  attain  higher  signal  to  noise  ratio.  Figure  2.2  illustrates  this  concept  where  a 
set  of  K  measurements  {g2}  are  taken  with  a  corresponding  K  set  of  optimal  aperture  codes  {w2} 
where  i  indexes  the  ith  measurement.  The  mapping  B  reorders  the  set  of  measurements  {g2}  into  the 
V  <  K  sets  {g1}  containing  only  information  of  the  spectral  bands  given  by  A^.  Each  subset  { g 2}  is 
sufficient  to  independently  reconstruct  the  desired  spectral  band  subset  {Fa,}-  The  new  coded  aperture 
agile  spectral  imaging  (CAASI)  system  described  in  this  Chapter  allows  multiple  measurement  shots 
and,  consequently,  it  has  the  aforementioned  advantages  over  CASSI.  The  new  capabilities,  however, 
come  at  the  expense  of  more  complexity  in  the  hardware  implementation.  Multishot  measurements  can 
be  attained  by  successively  shifting,  along  the  horizontal  axis,  the  fixed  coded  aperture  in  CASSI.  A 
novel  piezo-electrical  implementation  of  this  approach  is  described  in  [?,  ?].  An  approach  more  suitable 
for  dynamic  selectivity  in  real-time  is  to  implement  the  time-varying  aperture  codes  on  a  spatial  light 
modulator  [?,  ?,  ?,  ?].  Device  implementation,  calibration,  wavelength  range  of  operation,  and  other 
hardware  considerations  are  critical  but  fall  outside  the  scope  of  this  Chapter.  These  will  be  considered 
in  depth  in  a  separate  publication.  The  methods  developed  here  are  mathematical  in  nature  and  thus  are 
applicable  on  any  hardware  platform  capable  of  acquiring  multiple  shots,  each  with  a  distinct  aperture 
code  pattern. 

2.2  Coded  Aperture  Characterization  in  CASSI 

The  coded  aperture  single  shot  spectral  imaging  system  is  depicted  in  Fig.  2.3  [?].  The  coding  is  realized 
by  the  coded  aperture  T(x,y)  as  applied  to  the  spatio-spectral  density  source  fo(x,y]\)  where  (x,y) 
are  the  spatial  coordinates  and  A  is  the  wavelength  resulting  in  the  coded  field  fi(x,y,\).  The  coded 
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Figure  2.3:  Representation  of  the  CASSI. 


density  is  spectrally  dispersed  by  a  dispersive  element  before  it  impinges  on  the  focal  plane  array  (FPA) 
as  f2{x,y,  A), 

f2(x,  y,X)=  j  J T(x',  y')f0(x',  y' ,  A )h{x'  -  aX  -  x,  y'  -  y)dx'dy'  (2.1) 

where  T(x',yf)  is  the  transmission  function  representing  the  coded  aperture,  h(xf  —  aX  —  x,y'  —  y)  is 
the  optical  impulse  response  of  the  system,  and  aX  is  the  dispersion  induced  by  the  prism  assuming  a 
linear  dispersion.  Each  spectral  slice  of  the  data  cube  fo(x,  y ;  A)  is  thus  spatially  modulated  by  the  coded 
aperture  and  dispersed  by  the  dispersive  element  [?].  The  compressive  measurements  across  the  FPA  are 
realized  by  the  integration  of  the  field  /^(x,  y,  A)  over  the  detector’s  spectral  range  sensitivity.  The  source 
/o(V,?/,A)  can  be  written  in  discrete  form  as  (F/C)nm  where  n  and  m  index  the  spatial  coordinates, 
and  k  determines  the  kth  spectral  plane.  Assuming  that  the  band  pass  filter  of  the  instrument  limits  the 
spectral  components  between  Ai  and  A2  and  the  side  length  of  the  square  detector  pixel  is  A the  number 
of  resolvable  bands  L  is  limited  by  L  —  a*2^1 .  The  spectral  resolution  is  limited  by  Additionally, 
it  is  assumed  that  the  side  length  of  the  coded  aperture  square  pixel  Ac  satisfies  kAc  =  A where  k  >  1 
is  an  integer.  The  horizontal  and  vertical  spatial  resolutions  are  thus  limited  by  A The  height  N  of 
the  resolvable  data  cube  is  limited  by  N  =  min(^S  where  Uh  and  Vh  are  the  vertical  physical 
dimensions  of  the  coded  aperture  and  the  detector  respectively.  The  width  M  of  the  resolvable  data  cube 
is  limited  by  M  =  min(^^,  ^  —  L  +  1)  where  Uw  and  Vw  are  the  horizontal  physical  dimensions  of 
the  coded  aperture  and  the  detector  respectively.  Following  the  mathematical  model  in  [?],  the  coding  is 
realized  by  an  aperture  pattern  (T)nm. 


The  compressed  sensing  measurements  at  the  focal  plane  array  can  be  written  in  the  following  matrix 
form: 

L 

(^)nm  ^  ^(Ffc)n?m+fc (T)n;m_|_fe  T  (CU  )n,m  (2.2) 

k=l 


where  u  represents  white  noise,  L  is  the  number  of  spectral  bands,  and  n,  m  index  the  pixels  on  the 
detector.  A  typical  example  of  the  measurement  process  using  (2.2)  is  shown  in  Fig.  2.1.  The  expression 
in  (2.2)  can  be  expressed  as 

g  =  Hf  =  +  (2.3) 


where  g  and  f  are  vector  representations  of  G  and  F,  respectively  [?].  These  vector  representations 
will  be  described  shortly.  H  is  the  projection  matrix  that  takes  into  account  the  effects  of  the  coded 
aperture  and  the  dispersive  element.  T  is  a  Kronecker  basis  representation,  and  6  is  a  sparse  coefficient 
vector  representing  f.  The  characteristics  of  H  in  relation  to  incoherence  and  the  restricted  isometry 
property  (RIP),  needed  for  signal  reconstruction  in  a  set  of  undetermined  system  of  linear  equations  are 
studied  in  [?].  If  the  aperture  code  pattern  is  fixed  and  only  one  snap-shot  is  detected,  the  resultant 
spectral  imager  is  the  so-called  single  disperser,  or  single-shot,  CASSI  architecture  [?].  In  this  case, 
the  entire  3-D  multispectral  image  cube  is  compressed  into  a  single  2-D  compressive  image  measure¬ 
ment  at  the  FPA.  The  spectral  image  cube  f  can  be  reconstructed  by  solving  the  optimization  problem 
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Figure  2.4:  Model  equations  for  a  row  measurement.  A  slice  of  the  data  cube  fg  impinges  on  a  row  of 
the  coded  aperture  tl  =  r  o  w2  to  produce  the  coded  slice  The  elements  of  the  coded  slice  are 

reordered  into  the  matrix  Tl.  The  detector  output  g*  =  T1ul  sums  the  columns  of  TL 


f  =  Tjargmin^/  ||g  —  +  t||0,||i}  where  r  >  0  is  a  regularization  parameter  that  balances  the  con¬ 

flicting  tasks  of  minimizing  the  least  square  of  the  residuals  while,  at  the  same  time,  yielding  a  sparse 
solution.  Several  algorithms  exist  in  the  compressive  sensing  literature  to  solve  this  inverse  problem 
[?,  ?,  ?,  ?]. 

H  in  (2.3)  is  an  N  ■  (M  +  L  —  1)  x  N  •  M  •  L  matrix  mapping  the  vector  f ,  which  indexes  the  3D 
data  cube,  into  the  vector  g.  Note  from  (2.2)  and  (2.3)  that  the  matrix  H  embeds  both,  the  coded 
aperture  T  and  the  shifting  operation  of  the  dispersive  element.  Since  the  aperture  code  is  hidden  in 
H  and  since  it  is  not  directly  observed  in  (2.2)  and  (2.3),  its  optimization  is  difficult  using  this  CASSI 
formulation.  Given  that  this  analysis  only  considers  the  dispersive  and  the  coded  aperture  effects,  the 
matrix  H  is  binary;  however,  due  to  non-linearity  in  the  dispersive  element,  non  ideal  optical  instruments 
and  misalignments  between  the  detector  pixels  and  the  spatial  light  modulator,  these  matrices  have  real 
values  in  practice  [?].  For  the  purposes  of  this  Chapter  a  binary  H  is  used.  Additionally,  the  model  used 
does  not  consider  noise.  The  modifications  on  H  to  address  the  physical  limitations  of  the  optical  system 
can  be  readily  done  but  omitted  here  in  order  to  better  illustrate  the  aperture  code  optimization  problem. 

Since  the  dispersive  element  shifts  the  propagating  light  along  the  horizontal  dimension  only,  the 
spectral  coding  in  CASSI  is  independent  from  one  row  to  another.  Hence,  it  is  possible  to  construct  an 
aperture  code  model  for  just  one  row  of  the  measurement  matrix,  which  is  then  replicated  for  all  rows  of 
the  detector.  To  this  end,  define  the  vector  fg  that  indexes  one  slice  of  the  data  cube  (F k)nm  when  the 
index  n  is  fixed  to  q.  More  specifically, 

(fg)(/c— i)-M+m  =  (Ffc)gm,  for  m  =  1, . . . ,  M  fc  =  l,...,L.  (2.4) 

Note  that  the  elements  of  fg  are  the  only  elements  of  the  source  f  having  an  impact  on  the  qth  row  of  G. 
Let  the  qth  row  of  G  be  gr  The  projection  of  fq  into  gq  is  then  represented  as 

g*  =  Hgrfgr.  (2.5) 

for  q  =  1, _ ,7V.  Figure  2.4  illustrates  the  various  vector  representations  in  (2.4)  and  (2.5)  of  the 

source  data  cube.  Since  the  spectral  coding  along  rows  of  G  are  mutually  independent,  the  vectors  gg, 
can  be  arranged  as  in  (2.3)  where  g  and  H  are  unfolded  showing  the  mutual  independence  of  the  row 
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where  the  matrices  0  have  (M  +  L  —  1)  x  M  •  L  zero  elements,  the  NM L-long  vector  f  =  [f^  fj  ... 


is  the  concatenation  of  the  vectors  fq  in  (2.4),  and  the  vector  g  =  [gi  g2  •  •  •  gw]T  is  the  concatenation 
of  the  row  vectors  gq  representing  the  qth  row  at  the  detector.  The  matrix  Hq  represents  the  effect  of 
the  dispersive  element,  the  coded  aperture,  and  the  integration  on  the  detector  for  the  qth  row  of  the 
measurement  matrix  G.  The  Hq  matrix  is  given  by 


(2.7) 


where  diag{(T)q?i,  (T)qj2,  • . . ,  (T)q?M}  is  an  M  x  M  diagonal  matrix  with  the  elements  of  the  qth  row  of  T 
in  its  diagonal.  Observe  that  the  specific  structure  of  Hq  determines  the  structure  of  H  in  (2.6)  and  (2.3). 
In  order  to  separate  the  effects  of  the  various  optical  elements,  the  representation  Hq  can  be  expressed  as 
H q  =  D T,  where  the  matrix  D  represents  the  mapping  operation  of  the  dispersive  element.  The  matrix 
D  has  at  most  L  non  zero  elements  in  each  row.  The  specific  structure  of  the  (M  +  L  —  1)  x  ML  matrix 
D  is  given  by 


diag{(T)qji, . . . ,  (T)q  M} 

OlxM 

OlxM 

OlxM 

diag{(T)q>i,...,(T)q>M} 

OlxM 

OlxM 

OlxM 

OlxM 

diag{(T)q,1} . . . ,  (T)q,M> 

D  = 


I M 

OlxM 

OlxM 

OlxM 

I M 

OlxM 

OlxM 

OlxM 

I M 

(2.8) 


where  1m  is  an  M  x  M  identity  matrix  and  Oi  xm  are  zero- valued  row  vectors.  The  matrix  T  represents 
the  effect  of  the  coded  aperture  T  on  fq  and  has  the  following  structure 


T  =  II  ®  diag{(T)q;i,  (T)q?2,  •  •  • ,  (T)q;M} 


(2.9) 


where  II  is  an  L  x  L  identity  matrix.  T  is  thus  an  M  •  L  x  M  •  L  diagonal  matrix  whose  elements  are  L 
repetitions  of  the  elements  in  the  qth  row  of  T.  The  mapping  of  fq  into  gq  in  (2.5)  is  then  written  as 


g  q  —  q 


(2.10) 


The  vector  representation  in  (2.10),  characterizing  CASSI,  will  be  critical  in  formulating  the  aperture 
code  optimization  problem  described  in  the  following  sections. 


2.3  Coded  Aperture  Agile  Spectral  Imaging  System  (CAASI) 

The  spectral  selectivity  is  determined  by  the  aperture  codes  used  in  each  shot.  The  mathematical  model 
of  the  ith  measurement  shot  of  the  CAASI  system  is  identical  to  that  shown  in  (2.2)  for  the  CASSI  system 

L 

(G*)n?m  =  ^^(F/c)n?m+/c(T2)n?m+/c  +  ( Wl)n,m  (2.11) 

k= 1 

for  i  =  {1, . . . ,  AT},  where  K  is  the  number  of  shots,  except  that  the  ith  aperture  code  pattern  used 
to  sense  G2  is  different  from  one  shot  to  another.  Notice  that  the  coded  aperture  Tl  can  be  seen  as 
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two  separate  coded  apertures  applied  in  tandem.  The  first  coded  aperture  R  is  a  random  binary  code 
necessary  to  attain  randomized  measurement  in  CASSI.  The  second  coded  aperture  W2  is  the  ith  code 
optimized  to  achieve  a  specific  spectral  band  selectivity.  T2  is  thus  represented  as  the  Hadamard  product 
(T2)n?m  =  (R)n?m(W2)n?m  where  W2  is  the  time  varying  element  in  each  CAASI  measurement.  In  matrix 
form,  the  CAASI  system  of  equations  results  in 

g2  =  H2f  +  ui  i  G  1  ,...,AT.  (2.12) 

As  in  the  single  shot  CASSI  system,  the  row  measurements  g2  are  decoupled  from  one  another  such 
that  the  H2  matrices  in  (2.12)  are  block  diagonal.  The  qth  row  of  G2  is  then  given  by  g2  =  H2f g,  or 
equivalently 

s'q  =  Dff,.  (2.13) 

Note  in  (2.13)  that  the  source  vector  ^q  does  not  contain  the  index  i,  since  we  assume  that  the  source 
remains  stationary  during  the  time  interval  when  the  K  shots  are  measured.  The  only  difference  at  this 
point,  between  the  static  case  of  the  CASSI  equation  in  (2.10)  and  the  CAASI  equation  in  (2.13),  is  that 
the  matrix  T2  changes  with  each  shot.  The  matrices  T2  have  the  structure  given  in  (2.9),  but  in  this 
case,  the  qth  row  of  T2  can  be  expressed  as  t2  =  r  o  w2  where  o  denotes  the  Hadamard  multiplication  of 
r  and  w2,  the  qth  row  of  the  matrices  R  and  W2,  respectively.  Figure  2.4  depicts  the  interaction  of  the 
vector  t2  with  the  incoming  wave  fq.  In  the  static  case  of  (2.9),  T  does  not  have  a  time  varying  term  and 
thus  t  =  r.  The  matrix  T2  can  be  expressed  as  T2  =  II  0  diag{(w2)i(r)i, . . . ,  (w2) m (r) m} •  Using  the 
product  property  of  the  Kronecker  product,  T2  is  written  as 


Tl  =  (II  <g>  diag{(w*)i, . . . ,  (w')m})  (\l  <g>  diag{(r)i, . . . ,  (r)M})  (2.14) 

^ - * -  S 

W2  ti 

where  W2  represents  the  component  of  the  coded  aperture  w2  and  7 Zl  reflects  the  effect  of  r.  While 
the  random  matrix  7 Z  is  not  tunable  in  this  design,  the  matrix  W2  contains  weight  elements  of  the 
coded  aperture  that  will  be  optimized.  H  is  attained  via  Bernoulli  random  realizations,  however,  random 
matrices  such  as  a  scrambled  Hadamard  can  be  used  to  construct  this  random  component.  Using  (2.14) 
in  (2.13)  we  obtain  g2  =  DW27Jfg.  Letting  p  =  VAq  be  the  effect  of  the  vector  rq  on  the  source  fg,  then 
g2  can  be  written  as 

si  =  D  Wp.  (2.15) 

Making  some  matrix  indexes  manipulations,  Eq.  (2.15)  can  be  written  as 


where  ul  is  the  L  long  vector  wl 


g*  =r2uL 

and  the  matrix  T2  is  given  by 


Pl(wl)i  P2(wz)2  ... 

0  Pm+i  (wz)i  . .  . 

0  0 

0  0  0 


Pm+l-i(w*)L-i. 

P2M-\-L  —  2  (w*)l-2 
P(L-1)M  +  1  (wTl 


Pm(w1)m 

P2M-i(w')m_i 
P3M-2  (wZ)m-2 


0 

P2M  (w7, )  m 

P3M-i(w2)m-i 


pLM-(L-l)  (WZ) 


M-L  +  l 


(2.16) 


o 

o 

Plm(w*)m 

(2.17) 


where  pj  is  the  jth  element  of  p.  It  can  be  observed  that  the  (w l)i  terms  have  a  circular  shift  on  l  within 
each  row  of  T\  Additionally,  only  L  pixels  in  the  detector  are  affected  by  a  point  source  in  the  data  cube. 
Hence,  only  L  elements  of  w2  in  (2.17)  need  to  be  optimized.  Define  the  L-long  sub-segments  of  w2  as 
w2,  such  that  w2  is  formed  by  the  modulo  repetition  operator 


(w2)m  =  (w l)mod(m,L)  for  77i  =  1, . . . ,  M  and  L  <  M. 

Equation  (2.18)  can  be  expressed  in  matrix  form  as 

W2  =  U M'  0  W2 


(2.18) 


(2.19) 
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where  M'  =  ^  is  the  length  of  the  one  element  vector  u m'-  If  the  integral  part  of  ^  is  not  an  integer 
then  the  first  elements  of  w2  are  given  by  (2.19)  and  the  remaining  by,  (w2)m  =  (w l)m-L,  for 

m  =  +  1, . . . ,  M.  Reordering  the  elements  of  T2  in  (2.17)  such  that  the  (w l)i  terms  appear  in 

ascending  order,  and  using  the  periodicity  of  w2  given  in  (2.18),  we  obtain 


r  = 


Pl(w*)i  pM+l(wl)i  ...  P(L-l)M  +  l(wl)i  •••  Plm-(l- l)(W")l  0 

0  P2(w*)2  •••  P(i-2)M  +  2(w2)2  •  •  •  P(L-l)M-(L-2)  (w2)2  pLM-{L-2)  (wl)2 


n  T 


0 


0 


’•  P(I>-3)M  +  3(W2)3  •  •  •  P(JL-2)M-(I,-3)  (WZ)3  P(L  —  1) M  —  (L  —3)  (w2 ) 3 


Pm(w1)l 


P2M  (w* )  _L 


...  Plm{w1)l. 


.  (2.20) 


Note  that  the  matrix  product  r2uz,  in  (2.16)  and  the  product  r2uz,  of  the  reordered  matrices  in  (2.20) 
are  identical.  Hence  (2.16)  can  be  rewritten  as  g2  =  T1ul.  The  expression  for  g2  can  be  further  simplified 
by  factoring  out  the  elements  w2  in  (2.20)  leading  to 


si  =  X  w1  (2.21) 

where 


Pi 

Pm+ 1 

P(L-l)Af+l 

PLM-(L-l) 

0 

0 

0 

P2 

P(L-2)M+2 

P(L-l)M-(L-2) 

PLM—(L—2) 

0 

X  = 

0 

0 

P(L—3)M+3 

P(L  —  2)M—(L—3) 

P(L  —  l)M—(L—3) 

0 

0 

0 

PL 

Pm 

P2M 

Pl{m) 

is  a  matrix  containing  the  spectral  source  elements  multiplied  by  the  random  component  of  the  aperture 
code,  which  are  then  spectrally  shifted.  Equation  (2.21)  is  a  compact  expression  for  g2 ,  where  w2  describes 
the  effect  of  the  optimizable  part  of  the  coded  aperture.  The  effect  of  the  others  optical  elements  on  the 
data  cube  is  accounted  for  in  the  reordering  of  the  data  matrix  X.  The  complete  vector  g2  is  now 
succinctly  expressed  as 


■  si 

■  xxwi  ■ 

.  giv  . 

Xjvw^ 

where  the  structure  of  the  matrices  Xj  accounting  for  the  jth  data  cube  slice  is  given  by  (2.22),  and 
w2  is  the  designable  component  of  the  qth  row  of  the  coded  aperture.  In  matrix  notation  (2.23)  can  be 
expressed  as 

g2  =  XW  (2.24) 

where  W*  =  {w^T,  w^T  . . .  and  X  =  diag{Xi,  X2, . . . ,  X^v}  is  a  block  diagonal  matrix.  Equa¬ 

tions  (2.21)  and  (2.24)  provide  a  new  representation  for  the  measurements  of  CAASI  when  multiple  shots 
are  admitted.  The  new  formulation  is  illustrated  in  Fig.  2.5.  In  the  traditional  model  (Fig.  2.5(a)),  the 
input  is  the  data  cube  f  and  the  output  g2  is  characterized  by  the  matrix  H2.  In  the  new  representation 
(Fig.  2.5(b)),  the  input  is  the  matrix  X  and  the  output  is  obtained  by  varying  the  matrix  W\  The 
matrix  X  represents  a  reordered  version  of  the  vector  f  so  as  to  account  for  the  spectral  sheering  and 
R,  the  random  term  of  the  coded  aperture.  Note  that  the  elements  of  W*  in  the  new  representation  in 
(2.24)  explicitly  contains  the  components  of  the  coded  apertures  that  must  be  optimized.  Observe  that 
the  matrices  H2  in  (2.12)  have  a  total  of  ( N)(M  +  L  —  1)  x  (TV  •  M  •  L)  elements  each,  compared  to  the 
N(M  +  L  —  1)  x  L  elements  of  the  matrix  X  in  (2.23).  The  new  model  is  thus  significantly  more  compact 
than  that  of  (2.12)  and  consequently  it  is  better  suited  for  iterative  optimization  algorithms. 
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Figure  2.5:  (a)  The  traditional  model  where  the  data  cube  f  is  processed  with  the  highly  sparse  matrix 
H2.  The  coded  aperture  pattern  is  hidden  in  H2;  (b)  new  model  of  CAASI:  f  is  first  re-ordered  and 
expanded  into  the  matrix  X  which  is  then  processed  by  the  weight  matrix  W2  whose  elements  are  the 
coded  aperture  patterns. 

2.4  Aperture  Code  Optimization  in  Multishot  CAASI 

Each  measurement  in  CAASI  uses  an  optimized  W2  coded  aperture  component.  Since  the  spectral  coding 
is  mutually  independent  between  the  rows  of  the  detector,  it  is  sufficient  to  design  one  row  of  W2.  The 
remaining  rows  are  completed  by  shifted  replicas  of  the  first  such  that  W2  =  uatC)w2.  Furthermore,  since 
w2  is  a  concatenation  of  several  identical  segments  w2,  the  optimization  of  the  K  W2  coded  aperture 
matrices  reduces  to  the  optimization  of  K  w2  vectors. 

Before  the  optimization  problem  is  formulated,  it  is  useful  to  first  simplify  the  notation  by  deleting 
the  subindex  q  in  the  above  formulation,  which  indicates  a  particular  slice  of  the  CASSI  system.  It  is 
also  useful  to  note  the  similarities  of  the  CAASI  system  to  a  multichannel  filtering  problem  where  the 
input  signal  is  the  data  cube  f  and  the  filter  coefficients  are  the  elements  of  the  coded  aperture  w2.  The 
filtering  problem,  in  this  case,  aims  at  optimizing  the  coefficients  (coded  apertures)  such  that  the  linear 
combination  of  the  outputs  at  each  jth  position  at  the  detector  is  as  close  as  possible  to  a  desired  response, 
one  that  only  contains  information  of  the  specific  bands  of  interest.  The  coded  aperture  design  problem 
is  thus  analogous  to  finding  the  set  of  the  optimal  filter  coefficients  in  filter  design.  More  formally,  given 
an  input  signal  fs  and  a  desired  output  signal  d,  then  the  objective  is  to  find  the  K  optimal  filter  weight 
vectors  (aperture  codes  w1,  w2,...,wK)  and  a  set  of  associated  weights  bij  such  that  the  error  vector 
e  =  [ei,e2, . . . ,  cm+l- i]  is  minimized  according  to  an  error  criterion.  Figure  2.6  depicts  the  aperture 
code  optimization  problem  as  a  filtering  problem  for  the  jth  position  on  a  detector  row.  To  proceed  with 
the  filter  design  approach,  a  test  image  cube  fs  is  first  created  as  follows.  Define  a  sequence  of  distinct 
prime  numbers  {si,  52, . . . ,  sl}  such  that  si  ^  Sj  for  all  l  ^  j.  Choosing  the  set  of  primes  properly,  it  is 
possible  to  guarantee  that 


L 

^  ^  sj zj  7 ^  si  (2.25) 

3  = 1 
3& 

for  all  /,  where  the  coefficients  Zj  are  binary,  Zj  G  {0,1}.  If  a  coefficient  set  {zj}j '=1  is  different  from 
another  in  at  least  one  element,  then  from  the  properties  of  prime  numbers  it  follows  that 

Y^j=izjsj  7 ^  Xq=i  z'jsr  Given  the  set  of  prime  numbers,  the  test  data  cube  to  be  used  in  the  filter 
design  is  then  set  to 

(F^)nm  =  sk  for  all  m,  n.  (2.26) 

Hence  Fs,  is  a  cube  where  each  spectral  slice  takes  on  the  constant  value  of  sk  across  all  the  spatial 
positions  m  and  n.  A  slice  of  the  data  cube  at  n  =  q  in  vector  form  is  denoted  by  fs  across  all  m  and  k 
.  Clearly  (F|  )n?m  ^  (F^2)n?m  if  k\  ^  The  mapping  operation  converting  fs  into  X  is  shown  in  Fig. 
2.6  as  the  re-indexing  operation. 
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Figure  2.6:  Code  aperture  optimization  as  a  filter  design  problem.  The  input  is  a  row  of  the  data  cube 
fs  and  the  desired  signal  is  f°.  The  filter  coefficients  (aperture  codes  w^w2,. .  .,wx)  and  the  coefficients 
bij  are  optimized  by  minimizing  a  cost  function  of  the  errors  ei,  e2, . . . ,  6m+l- i-  The  error  at  given  jth 
position  at  a  row  at  the  detector  is  the  difference  between  the  linear  combination  of  the  measurements 
{Q)j  and  (d)j. 


From  (2.25)  and  the  CAASI  equation  (2.21),  given  the  detector  measurements  sequence  at  the  detector, 
g2,  it  is  possible  to  determine  which  spectral  bands  are  present  in  the  linear  combination  yielding  g2. 
This  is  the  case  since  (F%)qm  and  X  are  known.  The  desired  signal  d,  shown  in  Fig.  2.6,  is  created  from 
Fd  which  contains  only  the  desired  bands  and  is  given  by 


for  all  m,  n  and  where  the  L-long  binary  vector  A  indexes  the  subset  of  spectral  channels  of  interest.  In 
particular  (A)&  is  defined  as 


Wk 


1  if  the  kth  spectral  band  is  of  interest 
0  otherwise, 


(2.28) 


for  k  G  {1 The  qth  slice  of  Fd  is  fd  which  is  thus  forced  to  have  only  the  spectral  components  of 
interest.  The  output  at  the  detector,  when  the  desired  data  cube  fd  is  sensed,  is  given  by 

d*  =  Xdw2  (2.29) 

where  Xd  is  a  matrix  whose  elements  are  obtained  by  the  re-indexing  of  fd.  Since  Xd  contains  only 
information  from  the  desired  spectral  bands,  then  the  desired  signal  d2  can  be  calculated  by  a  code  w2 
that  simply  adds  the  row  elements  of  Xd.  This  operation  is  analogous  to  using  an  “all  pass”  code  w2 
that  permits  to  pass  all  the  spectral  information  from  Xd.  In  this  case,  the  all-pass  code  w2  =  u^,  in 
effect,  simply  sums  all  elements  of  Xd.  The  vector  d2  in  (2.29)  does  not  depend  of  the  coded  aperture  w2 
and  it  is  the  desired  output  d.  Note  that  the  desired  signal  d  is  unique  in  the  sense  that  only  the  specific 
combination  of  spectral  bands  given  by  A  can  generate  the  output  of  CAASI  equal  to  d.  This  is  precisely 
the  motivation  of  using  the  set  of  prime  numbers  in  (2.25)  as  elements  of  the  data  cube  fs.  To  further 
illustrate  the  underlying  concepts,  Fig.  2.7  shows  a  typical  slice  of  the  data  cube  fs  and  the  corresponding 
structure  of  the  matrix  X.  Observe  that  the  elements  of  the  output  g2  are  calculated  by  the  Hadamard 
product  of  the  coded  aperture  w2  and  each  row  of  X.  In  Fig.  2.6  the  operation  of  extracting  Xj  the  jth 
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Figure  2.7:  A  slice  of  the  data  cube  fs  is  modulated  by  the  vector  r  and  then  reorganized  into  the  matrix 
X.  The  output  at  the  detector  is  calculated  as  g2  =  Xw\ 


row  of  X  is  indicated  with  the  “Row-Select”  operation.  From  Fig.  2.6,  the  error  at  jth  position  in  the 
detector  is 

ej  =  (d )j  ~  0 Q%  j  =  1, . . . ,  M  +  L  -  1  (2.30) 

where  (g)j  =  J2f=i  bijX-J w\  Notice  that  d  and  Xj  are  known  terms  and  w2  represents  the  unknown  aper¬ 
ture  codes  pattern  to  be  designed.  The  coefficients  bij  permit  to  use  all  the  information  available  from 
the  K  shots  to  minimize  the  jth  component  at  a  row  at  the  detector.  The  coefficients  bij  are  restricted 
to  the  three  values  {  —  1,0, 1}  that  the  value  —1  indicates  the  subtraction  of  a  set  of  bands,  0  indicates 
the  suppression  of  a  set  of  bands,  and  1  indicates  the  summation  of  a  set  of  bands.  The  errors  origi¬ 
nated  from  the  measurements  in  (2.30)  are  concatenated  into  the  common  array  e'  =  [ei,  e2, . . . ,  cm+l-i\- 

The  component- wise  squared  error  metric  is  first  calculated  such  as  e  =  e'  o  e'.  The  jth  element  of 
e  can  be  expressed  as  a  function  of  the  coded  aperture  as  (e)j  =  ((d) j  —  Ylf=  i  .  Note  that 

(e)j  =  0  occurs  only  when  the  coded  apertures  are  able  to  exactly  extract  the  desired  bands  A  in  the  jth 
position  in  the  detector.  Note  also  that  when  (e)j  >  0  then  there  are  additional  spectral  bands  included, 
or  not  all  desired  bands  are  present  in  the  jth  entry  of  the  q  vector.  Given  that  the  spectral  selectivity 
consists  on  obtaining  at  each  position  in  the  detector  exclusively  the  desired  spectral  components,  it  is 
equivalent  to  maximizing  the  number  of  elements  (e)j  equal  to  zero.  In  other  words,  it  is  equivalent  to 
minimizing  the  Lq  norm  of  the  vector  e.  Hence,  one  seeks  to  obtain  the  vectors  W  =  [w1, . . . ,  wK]  and 
the  coefficients  B  =  [{bn, . .  .,bK1}T, {&i(m+l-i),  •  •  • ,  bK(M+L-i)}T]  such  that 

argmin  ||e||0  (2.31) 

W.B 

subject  to  (W )ki  G  {0, 1}  k  =  1, . . . ,  L,  i  =  1, . . . ,  K 

(B )ij  e  {-1,0, 1}  j  =  l,...,M  +  L-l. 

Since  the  measurements  gz  contain  complementary  information  between  them,  the  functionality  of  the 
coefficients  bij  is  to  increase  the  number  of  solutions  in  (2.31)  such  that  an  iterative  algorithm  can 
converge  more  easily.  Since  the  elements  of  and  bij  are  binary  and  ternary,  respectively,  solving 
(2.31)  is  an  NP  hard  problem.  An  approximate  solution,  however,  can  be  sought  using  a  heuristic 
approximation  such  as  a  genetic  algorithm  (GA)  [?].  GAs  are  specially  well  suited  to  solve  optimization 
problems  that  involves  binary  variables.  This  optimization  technique  has  already  been  used  in  several 
optics  optimization  problems  [?,  ?].  The  genetic  algorithm  (GA)  approach  to  solve  (2.31)  is  shown  in 
Table  2.1.  The  stopping  criterion  to  solve  (2.31)  is  that  ||e||o  =  0.  Notice  that  this  stopping  criterion 
conduces  to  a  solution  where  the  optimal  codes  lead  to  measurements  containing  exclusively  the  desired 
bands.  A  measurement  containing  the  desired  components  can  be  constructed  from  the  results  of  the 
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optimization  as 


Q  = 


K 


K 


K 


i= 1 


i=  1 


i= 1 


(2.32) 


where  B*  contains  the  optimal  linear  coefficients  and  w*2  are  the  optimal  coded  apertures  attained  in 
(2.31). 


Table  2.1:  Iterative  Optimization  of  coded  apertures  {w*1, . . . ,  w*K}  and  the  optimal  matrix  B*. 


Inputs 

Desired  spectral  bands  A,  Number  of  shots  K ,  Prime  numbers  S&,  Population  size  S , 
Number  of  iterations  A/",  Mutation  percentage  pm. 

Initialize 

Data  Cube  Sk  — >  fs  — *  Xs.  Desired  Signal  (s*;,  A)  — >  f d  — >  Xd  — >  — yd  — >  T>  =  d  ®  u^. 

B  =  {b1 , . .  . ,  b^},  where  1C  =  3K ,  and  bmi  ^  bm2  for  each  m\  ^  m2,  m\  =  1, . . . ,  1C, 
and  m2  =  1, ...  ,JC. 

Global  minimum  £*  =  M  +  L  —  1.  Population:  For  l  —  1, .  . .  ,  S  compute: 

=  {wj , . . .  ,  wf  },  where  (w^)^.=Bernoulli(0.5). 

Iteration 

While,  n  <  N  and  £*  >  0 

For  £  =  1, .  . . ,  S  compute 

G(p  =  XP<n),  QP  =  G E<n)  =  (v  -  £)n))  o(v~ 

(4n)t  =  min  (E<n))  I  (£(n))<!  =  \\4n)\\0,(CM)ej  =  argmin  (e^)  . 

c  m  \  c  / jm  m  \  c  J jm 

End  For 

£*  =  arg  min^71))^. 

If  (€(n)h*  <e,  then  : 

T  =  «{n))/.,{W1,...Iw*^  =  P(”),  (B %•  =  (B)fc(£(»)W 

End  If 

{P^n+1), . . .  ,  P^n+1)}  =  Mu  (Cr  (SI  ({P^n),  •  •  • ,  P£°},£(n)))). 

End  While 

Output 

Optimal  aperture  codes  w*1, .  . .  ,w*K. 

Optimal  linear  coefficient  B*. 

Note 

j  =  l,...,M  +  L  —  1  and  k  =  1, . . . ,  L. 

A  property  of  the  optimized  aperture  codes  is  that  if  the  codes  w*1, . . . ,  w*K  and  their  respective 
optimal  coefficient  matrix  B*  can  extract  the  spectral  bands  given  by  A,  then  the  same  codes  with  a 
different  coefficient  matrix  Bn*  can  extract  the  spectral  bands  given  by  An  with 

(\n)k  =  W((k-n)L),  k  =  l,...L  (2.33) 

where  (k  —  ti)l  indicates  the  k  —  n  modulo  L  operation  or  equivalently  a  n  circular  shift  of  period  L. 
This  property  requires  that  A  has  some  circular  period  J\f  and  that  the  ratio  jj>  is  an  integer.  Finding 
the  new  matrix  Bn*  involves  simple  algebraic  operations.  The  optimized  aperture  codes  are  thus  circular 
shift  invariant  to  the  vector  A. 


2.5  Simulations  and  Results 

2.5.1  Coded  Aperture  Optimization 

To  test  the  methodology  developed  in  Section  2.4,  the  algorithm  in  Table  2.1  is  used  to  find  the  optimal 
aperture  codes  for  several  desired  spectral  band  sets  A,  for  increasing  number  of  measurement  shots.  The 
experiments  use  a  test  data  cube  Fs  with  L  =  24  spectral  bands  in  504nm-757nm  and  256  x  256  spatial 
dimensions.  A  set  of  24  prime  numbers  Sk  is  generated  such  that  (2.25)  is  satisfied.  The  other  parameters 
of  the  algorithm  in  Table  1  are:  population  size  S  =  100,  mutation  percentage  pm  =  0.01,  and  maximum 
number  of  iterations  A f  =  5000.  In  the  experiment,  the  following  desired  spectral  bands  A  are  used  in 
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the  optimization 


A1  =  [100010100000000000110000]  A2  =  [100010001000100010001000] 

A3  =  [101000101000101000101000]  A4  =  [110000001100000011000000] .  (2.34) 

From  numerous  simulations,  it  can  be  observed  that  the  algorithm  converges  for  K  <  L  when  the  vector 


Figure  2.8:  Performance  evolution  as  a  function  of  the  iterations  when  the  vectors  (a)  A2  and  (b)  A4 
in  (2.34)  are  used  as  input  of  the  GA. 

A  is  circularly  periodic.  When  A  is  aperiodic,  the  algorithm  converges  for  K  =  L.  Figure  2.8  shows  the 
performance  evolution  as  a  function  of  the  iterations  when  the  vectors  A2  and  A4  are  used  as  input 
of  the  AG.  Figure  2.9  illustrates  the  optimal  aperture  codes  w*2  when  the  vectors  A1,  A2,  A3,  and  A4 
are  used  as  input  to  the  optimization  algorithm.  Additionally,  Fig.  2.10  shows  an  128  x  128  matrix 
whose  elements  are  constructed  by  replicates  of  the  optimal  coded  apertures  and  coded  by  a  random 
component.  Each  row  of  these  matrices  is  constructed  by  selecting  randomly  an  optical  code  w\  such 
that  the  entire  matrix  contains  information  of  the  all  optimal  codes.  In  addition  to  the  optimal  codes, 
the  optimization  provides  the  optimal  coefficient  matrix  B*  used  to  construct  the  elements  of  the  optimal 
virtual  measurement  q* .  Figure  2.11  illustrates  the  optimal  matrix  B*  when  A4  is  used;  due  to  space 
limitations  only  a  portion  of  this  matrix  is  shown.  Note  that  the  patterns  in  Fig.  2.11  exhibit  periodicity, 
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except  in  the  beginning  and  end  which  exhibit  edge  structures  that  are  different  from  that  of  the  main 
region. 


Figure  2.10:  An  128  x  128  realization  of  the  optimal  coded  apertures  for  the  vectors  (a)  A1,  (b)  A2,  (c) 
A3,  (d)  A4  given  in  (2.34). 


2.5.2  Filter  Bank  Representation 

The  filter  bank  spectral  model  aims  at  reconstructing  the  entire  data  cube  from  the  V  compressive 
measurements  where  each  measurement  has  information  only  from  one  of  V  disjoint  spectral  subsets  with 
the  constraint  of  using  only  K  <  L  shots.  Each  subset  can  be  reconstructed  separately  and  the  results  can 
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Figure  2.11:  A  part  of  the  matrix  B*  is  shown.  The  Ith  column  of  B*  represents  the  optimal  coefficient  to 
construct  the  measurement  element  (g*)i.  The  vector  A4  is  used  as  input  of  the  optimization  algorithm. 


be  merged.  The  disjoint  spectral  subsets  are  specifically  designed  such  that  the  shift  invariance  property 
of  the  coded  apertures  developed  in  previous  sections  is  used.  More  specifically,  suppose  that  the  L  long 
vector  A  has  circular  period  AT,  define  U  as  the  maximum  number  of  consecutive  one  elements  in  a  period 
of  A,  and  define  V  as  V  =  \jp\  Then  the  complete  data  cube  can  be  reorganized  in  V  compressive 
measurements  with  information  only  from  the  bands  given  by  {Ao,  Ai, . . . ,  Ay-i}  where  the  elements  of 
each  subset  An  is  given  by 

(*»)j=  W((k-nU))L  k  =  l,...,L.  (2.35) 

The  filter  bank  representation  is  attained  as  follows:  (a)  generate  a  set  of  vectors  {Ao, . . . ,  Ay_i}  where 
each  vector  A^  is  a  shifted  version  of  Ao  =  A;  (b)  verify  that  Ao  V  Ai  V  . . .  V  Ay_i  =  ul  where  V  is  the 
boolean  sum;  (c)  The  spectral  components  given  by  the  different  A^  can  be  reconstructed  separately  and 
the  results  can  be  merged  with  the  other  reconstructions.  Since  [w*1, . . . ,  w*K]  are  the  optimal  codes 

and  B*  is  the  optimal  coefficient  matrix  to  recover  A,  then  each  combination  {[w*1, . . . ,  w*K]T  ,  B**} 
can  be  seen  as  a  compressive  filter  bank  to  recover  A^,  where  B**  is  a  reordered  version  of  the  matrix 
B*.  A  representation  of  this  technique  was  illustrated  in  Fig.  2.2.  Define  the  “modp”  A  pattern  which 
elements  are  given  by 

(A)fc  =  mod(fc,p)  fc  =  1, . . . ,  L  (2.36) 

where  mod  is  the  modulo  operation,  p  is  the  circular  period  of  A,  and  L  is  the  total  number  of  bands. 
Here,  it  is  assumed  that  the  ratio  ^  is  an  integer.  These  “modp”  patterns  can  be  used  to  decompose 
the  data  cube  in  the  subsets  {Ao,  Ai, . . .  , Ap_i}  where  Ao  is  equal  to  A  in  (2.36)  and  the  others  A^  are 
given  by  (2.35)  with  U  =  1.  Notice  that  the  A^  subsets  satisfy  the  above  (b)  condition  necessary  to  apply 
the  filter  bank  decomposition.  Decompose  the  data  cube  in  the  “modp”  subsets  requires  p  shots.  The 
respective  optimal  coded  apertures  to  select  the  modp  patterns  is  referred  to  as  “modp  filter  bank  coded 
apertures”.  For  example,  referring  to  Eq.  (2.34),  A2  is  a  mod4  pattern  and  the  codes  of  Fig.  2.9(b)  are 
a  mod4  filter  bank  coded  apertures. 

2.5.3  Reconstruction  Algorithm 

A  number  of  CS  reconstruction  algorithm  are  available  in  the  literature  [?,  ?,  ?,  ?].  The  Gradient 

Projection  for  Sparse  Reconstruction  algorithm  (GPSR)  [?]  was  used  to  reconstruct  the  spectral  data 

cubes  in  our  simulations  as  it  exhibits  faster  computational  speed.  This  algorithm  reconstructs  of  the  data 

cube  f  by  solving  the  minimization  f  =  T{argmin  ||g*  —  HT6>/||2  +  rl|0/||i}?  where  g*  are  the  optimized 

O' 

measurements.  The  operation  that  involves  the  sparse  matrix  H  are  implemented  through  the  indexes 
of  non-zero  elements  of  the  matrix  H  instead  of  using  it  directly.  The  base  representation  T  is  the 
Kronecker  product  of  three  basis  T  =  \EriG)\Er2(S)\k3,  where  the  combination  0  T2  is  the  2D- Wavelet 
Symmlet  8  basis  and  T3  is  the  Cosine  basis  [?].  A  data  cube  of  satellite  images  is  used,  it  has  24  spectral 
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channels  and  256  x  256  pixels  of  spatial  resolution.  Reconstruction  was  realized  using  a  desktop  with  3.06 
GHz  Intel  Core  2  Duo  processor,  and  4GByte  1067  MHz  DDR3  RAM.  The  simulation  tool  was  Matlab 
2007a  version  7.4.0  for  MacOS  version  10.6.4.  The  simulations  were  realized  using  double  numerical 
precision.  The  maximum  number  of  iterations  in  the  GPSR  algorithm  is  700  and  the  stopping  criterion 
as  recommend  in  [?]  is  the  objective  function.  Figure  2.12(b),  Fig.  2.12(c),  and  Fig.  2.12(d)  depict  the 
reconstruction  results  for  the  first  spectral  band  of  the  data  cube  applying  the  technique  of  compressive 
filter  bank.  The  mod2,  modi2,  and  mod24  filter  bank  coded  apertures  are  used  in  the  decomposition. 
Clearly  the  SNR  increases  with  the  number  of  shots  used. 


(a)  (b)PSNR  34.29  dB 


(c)  PSNR  42.4  dB  (d)  PSNR  50.7  dB 

Figure  2.12:  Reconstruction  of  the  first  spectral  band  of  the  24  spectral  band  data  cube,  (a)  Original; 
reconstruction  using  (b)  mod4  filter  bank  coded  apertures,  4  shots;  (c)  modi2  filter  bank  coded  apertures, 
12  shots;  (d)  mod24  filter  bank  coded  apertures,  24  shots. 

A  comparison  between  the  compressive  filter  bank  technique  and  a  random  multishot  approach  was 
realized.  The  random  multishot  technique  consists  on  using  random  coded  aperture  patterns  to  realize 
the  compressive  measurements  [?,?].  In  that  technique,  no  optimization  of  the  coded  apertures  was 
carried  out.  Figures  2.13(b),  2.13(c),  and  2.13(d)  show  the  reconstruction  of  the  same  first  band  of  the 
data  cube  using  the  random  codes  for  K  equal  to  2,  12  and  24  shots.  Figure  2.14  illustrates  the  general 
results  in  terms  of  SNR  for  the  multishot  and  filter  bank  approaches  as  a  function  of  the  number  of  shots 
K .  Notice  that  in  the  multishot  case,  the  complete  spectral  data  cube  is  reconstructed  at  once.  In  the 
filter  bank  approach,  on  the  other  hand,  the  complete  data  cube  is  divided  in  K  disjoints  subsets  where 
the  nth  subset  An  is  given  by 


(An)j 


1  for  j  =  n,  K  +  n,  2 K  +  n, . . . ,  (^  —  1  )K  +  n 
0  otherwise. 


(2.37) 


where  n  =  0, . . . ,  if  —  1.  Each  desired  subset  is  reconstructed  independently  from  others.  The  recon¬ 
structed  subset  bands  are  then  merged  and  the  total  SNR  is  calculated.  The  results  show  that  when 
the  number  of  shots  increases,  the  filter  bank  approach  is  superior  in  quality  of  reconstruction.  Figure 
2.15  shows  the  respective  reconstruction  time  for  the  same  coded  apertures  used  in  Fig.  2.14.  For  the 
reconstruction  technique  Random  Multishot,  the  reconstruction  time  decreases  with  increasing  number  of 
shots  because  when  more  shots  are  used  the  system  of  equations  are  less  undetermined  and  the  objective 
function  decreases  more  quickly  through  the  iterations.  In  the  filter  bank  case  the  reconstruction  times 
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(c)  PSNR  34.68  dB  (d)  PSNR  34.74  dB 


Figure  2.13:  Reconstruction  of  the  first  spectral  band  of  the  24  spectral  band  data  cube,  (a)  Original;  (b) 
random  coded  aperture,  4  shots;  (c)  random  coded  aperture,  12  shots;  and  (d)  random  coded  aperture, 
24  shots. 


Figure  2.14:  Mean  PSNR  for  the  reconstructed  data  cube  as  a  function  of  the  number  of  shots.  The 
techniques  of  random  multishot  and  compressive  modp  filter  bank  are  shown. 


are  shown  using  both  a  single  processor  and  a  processor  for  the  reconstruction  of  each  subset  of  bands. 
The  approach  using  several  processors  is  coined  “parallel  filter  bank” .  The  parallel  filter  bank  approach 
provides  up  to  two  orders  of  magnitude  faster  reconstruction  than  the  random  multishot  approach.  In 
the  second  experiment,  we  observe  the  spectral  selectivity  of  the  optimized  coded  apertures.  In  this  case 
the  reconstruction  of  the  1st  and  18st  bands  are  of  interest.  The  vector  A4  in  (2.34)  is  used  as  input  of  the 
optimization  algorithm.  In  this  experiment,  the  results  for  K  equal  to  4,  8,  12,  and  16  shots  are  shown  in 
Figs.  2.16(c)-2.16(d),  2.16(e)-2.16(f),  2.17(a)-2.17(b),  and  2.17(c)-2.17(d)  respectively.  Note  that  in  Fig. 
2.8(b)  shows  a  minimum  number  of  K  =  8  shots  for  convergence  of  the  optimization  algorithm.  However, 
Figs  2.16(c)-2.16(d)  show  that  it  is  possible  to  reconstruct  the  desired  spectral  bands  with  fewer  number 


43 


Figure  2.15:  Reconstruction  times  for  the  full  data  cube  as  a  function  of  the  number  of  shots  when  is  used 
(a)  pure  random  coded  apertures  (Random  Multishot);  (b)  modp  filter  bank  optimized  coded  apertures 
(Filter  Bank);  and  (c)  modp  filter  bank  coded  apertures  using  a  processor  for  each  subset  of  bands  (Filter 
Bank  Parallel). 


of  shots. 


2.6  Conclusions 

We  have  developed  the  underlying  theory  and  optimization  framework  for  a  new  Coded  Aperture  Agile 
Spectral  Imaging  System  (CAASI)  which  adds  the  spectral  selectivity  property  to  CASSI.  The  mathemat¬ 
ical  model  of  CAASI  expresses  the  output  signal  as  the  product  of  a  reordered  version  of  the  source  input 
and  the  coded  aperture.  Using  the  new  CAASI  representation,  an  optimization  algorithm  is  developed 
which  finds  the  optimal  coded  apertures  to  be  used  in  a  multishot  system  such  that  spectrally  selective 
compressive  measurements  are  attained.  The  optimization  aims  at  minimizing  the  Lq  norm  between  a 
data  cube  and  a  desired  cube  containing  only  the  spectral  bands  of  interest.  A  Genetic  Algorithm  solves 
the  Lq  optimization  problem  and  finds  the  minimum  number  of  shots  for  convergence.  The  optimization 
algorithm  converges  to  K  <  L  shots  when  the  subset  of  desired  bands  exhibits  some  circular  periodicity 
which  is  originated  because  any  assumption  is  made  over  the  random  component  of  the  coded  aperture. 
The  designed  coded  apertures  exhibits  a  circular  shift  invariance  property  that  permits  to  divide  the 
reconstruction  of  the  complete  data  cube  in  V  <  K  less  complex  reconstructions  problems.  This  property 
of  the  optimal  codes  can  be  seen  as  a  compressive  filter  bank  decomposition.  The  filter  bank  approach 
reconstructs  either  the  complete  data  cube,  or  a  selective  subset  of  bands,  with  higher  PSNR  and  takes 
up  to  two  orders  of  magnitude  less  processing  time  than  the  traditional  multishot  approach  with  random 
coded  apertures. 
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(e)  PSNR  37.30  dB  (f)  PSNR  36.57  dB 


Figure  2.16:  Reconstruction  of  the  1st  and  18th  spectral  band  of  the  24  band  data  cube,  (a)  Original  1st 
band  (b)  original  18th  band;  reconstruction  of  the  respective  band  using  the  vector  A4  in  (2.34).  (c),  (d) 
4  shots  ;  (e),  (f)  8  shots. 
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(c)  PSNR  39.59  dB  (d)  PSNR  39.63  dB 

Figure  2.17:  Reconstruction  of  the  1st  and  18th  spectral  bands  indicated  in  Fig.  2.16(a)  and  2.16(b). 
Reconstruction  of  the  respective  band  using  the  vector  A4  in  (2.34)  for  (c),  (d)  12  shots;  (e),  (f)  16  shots. 
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Chapter  3 


Rank  Minimization  Coded  Aperture 
Design  for  Spectrally  Selective 
Compressive  Imaging 

3.1  Introduction 

Consider  again  the  Code  Aperture  Snapshot  Spectral  Imaging  (CASSI)  system  that  allows  capturing 
spectral  imaging  information  of  a  3D  cube  with  just  a  single  2D  measurement  of  the  coded  and  spectrally 
dispersed  source  field  [?].  More  formally,  suppose  a  hyperspectral  signal  T  G  RiVxMxLj  or  its  vector 
representation  /  G  is  S  sparse  on  some  basis  such  that  /  =  sdP  can  be  approximated  by 

a  linear  combination  of  S  vectors  from  with  S  <C  (N  •  M  •  L).  Here,  N  x  M  represents  the  spatial 
dimensions  and  L  is  the  spectral  depth  of  the  image  cube.  Compressive  sensing  shows  that  /  can  be 
recovered  from  m  random  projections  with  high  probability  when  mn  >  S\og(N  •  M  •  L)  <C  (N  •  M  •  L). 
The  CASSI  projections  are  given  by  yt  =  H p,  where  H  =  is  an  N(M  +  L  —  1)  x  (TV  •  M  •  L)  matrix 

and  4>  is  a  random  measurement  matrix  determined  by  the  coded  apertures  and  the  dispersive  element 
used  in  CASSI. 

Recently,  the  CASSI  spectral  imaging  architecture  has  been  extended  to  admit  multiple  measurement 
shots  [?,  ?,  ?].  The  multiple  measurements  are  attained  as  separate  FPA  measurements,  each  with  a 
distinct  coded  aperture  that  remains  fixed  during  the  integration  time  of  the  detector.  There  are  several 
advantages  to  multiple  shots.  First,  the  number  of  compressive  measurements  in  CASSI  may  not  meet 
the  minimum  needed  for  adequate  reconstruction.  CS  dictates  that  the  number  of  measurements  must  be 
in  excess  of  5  log  (TV  •  M  •  L).  Failure  to  collect  a  sufficient  number  of  measurements  leads  to  inadequate 
signal  reconstruction.  With  each  FPA  shot,  CASSI  collects  N(M  +  L  —  1)  additional  measurements. 
For  spectrally  rich  scenes  or  very  detailed  spatial  scenes,  a  single  shot  CASSI  measurement  may  not 
provide  a  sufficient  number  of  compressive  measurements.  Increasing  the  number  of  measurement  shots 
will  multiply  the  number  of  measurements,  thus  rapidly  overcoming  such  limitations. 

A  second  advantage  to  multiple  shots  is  that  spectral  selectivity  can  be  attained  by  coded  aperture 
design.  Notably,  the  coded  aperture  patterns  can  be  designed  so  as  to  maximize  the  information  content 
on  a  pre- specified  subset  of  spectral  bands  of  particular  interest  in  the  CASSI  compressive  measurements. 
Spectral  selectivity  is  of  interest  in  many  applications,  including  wide-area  airborne  surveillance,  remote 
sensing,  and  tissue  spectroscopy  in  medicine.  The  optimal  spectral  bands  in  airborne  surveillance,  for 
instance,  depend  on  atmospheric  conditions,  time  of  day,  the  targets  of  interest,  and  the  background 
against  which  the  targets  are  viewed  [?].  In  these  applications,  the  spectral  signatures  of  interest  live  in 
a  spectral  band  subspace.  Efforts  placed  on  acquiring  the  entire  spectral  image  cube,  to  then  throw  away 
a  large  portion  of  this  data  is  wasteful  in  many  regards. 

Chapter  2  provided  the  first  approach  for  spectral  selectivity  in  multi-shot  CASSI;  however,  the  selec- 
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tive  spectral  profiles  were  limited  to  periodic  patterns  and  the  minimum  number  of  shots  was  restricted 
to  the  periodicity  of  the  spectral  pattern.  In  most  practical  applications,  however,  the  spectral  profiles 
of  interest  are  not  periodic  and  the  number  of  shots  is  restricted  by  the  application.  The  main  contri¬ 
bution  of  this  Chapter  is  the  development  of  a  more  general  and  more  effective  mathematical  framework 
for  multi-shot  CAS  SI  and  the  corresponding  algorithms  for  coded  aperture  optimization  that  allow  the 
reconstruction  of  arbitrary  subset  of  bands,  periodic  or  aperiodic,  whilst  minimizing  the  required  number 
of  shots. 

The  CAS  SI  model  is  first  developed  for  a  slice  at  the  detector.  It  is  then  used  to  build  a  generalized  2D 
CAS  SI  model.  The  output  at  a  slice  in  the  detector  is  expressed  as  the  product  of  a  matrix  accounting 
for  the  dispersed  spatio-spectral  source  density  and  a  vector  accounting  for  the  coded  aperture.  The 
coded  aperture  for  the  ith  measurement,  denoted  as  t2,  is  seen  as  two  separate  coded  apertures  applied  in 
tandem  which  simplifies  their  optimization.  The  first  coded  aperture  w2  is  a  structured  code  optimized 
to  attain  the  spectral  band  selectivity.  The  second  coded  aperture  r2  is  a  pseudorandom  binary  code 
necessary  to  attain  randomized  measurements  in  CASSI.  The  coded  aperture  used  in  each  measurement 
is  obtained  by  the  Hadamard  product  t2  =  w2  o  r2.  The  coded  aperture  optimization  is  divided  in  two 
parts.  First  the  tunable  components  w2  are  optimized  to  achieve  spectrally  selective  measurements.  The 
required  number  of  shots  K  is  dictated  by  the  rank  of  a  matrix  containing  the  w2  codes.  Second,  the 
pseudorandom  components  of  the  coded  apertures  r2  are  optimized  so  as  to  minimize  the  required  num¬ 
ber  of  shots.  This  is  achieved  by  minimizing  the  rank  of  the  matrix  containing  the  set  of  codes  t2.  The 
optimized  set  of  coded  apertures  leads  to  spectrally  selective  compressed  measurements  with  improved 
characteristics  for  meeting  Restricted  Isometry  Property  (RIP)  of  the  projections  used  in  CASSI.  The 
compressed  measurements  are  then  used  to  reconstruct  the  desired  spectral  cube  using  a  gradient  based 
reconstruction  algorithm  for  sparse  signals. 

Notation:  The  following  font  notation  is  used  hereafter.  Bold  uppercase  Roman  and  Greek  letters 
represent  matrices.  Bold  lower  case  Roman  and  Greek  letters  represent  column  vectors.  The  entries 
of  the  matrix  Y  are  or  (Y)^-  and  the  entries  of  the  vector  y  are  yj  or  (y)j.  Subindices,  upper 
indices  and  calligraphic  fonts  are  used  to  represent  distinct  variables.  Hence,  the  variables  y y  and  y y 
are  distinct  column  vectors.  The  variables  Y^  and  y^  for  different  i  are  distinct  matrices  and  column 
vectors  respectively.  yT  is  a  row  vector  and  YT  is  the  transpose  of  the  matrix  Y.  The  power  of  a  matrix 
is  presented  using  parenthesis,  thus  the  power  two  of  the  matrix  0^  is  (0l)2- 

3.2  Matrix-Based  CASSI  Modeling 

The  coded  aperture  single  shot  spectral  imaging  system  is  depicted  in  Fig.  2.3  [?].  The  coding  is  realized 
by  the  coded  aperture  T(x,y)  as  applied  to  the  spatio-spectral  density  source  fo(x,y,\)  where  (x,y) 
are  the  spatial  coordinates  and  A  is  the  wavelength  resulting  in  the  coded  field  /i(x,t/,  A).  The  coded 
density  is  spectrally  dispersed  by  a  dispersive  element  before  it  impinges  on  the  focal  plane  array  (FPA) 
as  f2(x,y,  A), 


h(x,y,  X)=J Jt(x' ,y')f0(x' ,y' ,  X)h(x'-aX-x,y' -y)dx'dy' 

where  T(xf,  y')  is  the  transmission  function  representing  the  coded  aperture,  h(x'  —  aX  —  x,  y'  —  y)  is  the 
optical  impulse  response  of  the  system,  and  aX  is  the  dispersion  induced  by  the  prism  assuming  a  linear 
dispersion.  Each  spectral  slice  of  the  data  cube  is  thus  spatially  modulated  by  the  coded  aperture  and 
dispersed  by  the  dispersive  element  [?].  The  compressive  measurements  across  the  FPA  are  realized  by 
the  integration  of  the  field  y ,  A)  over  the  detector’s  spectral  range  sensitivity.  The  source  fo(x ,  y ,  A) 
can  be  written  in  discrete  form  as  where  i  and  j  index  the  spatial  coordinates,  and  k  determines  the 
kth  spectral  plane.  Assuming  that  the  band-pass  filter  of  the  instrument  limits  the  spectral  components 
between  Ai  and  A2  and  the  side  length  of  the  square  detector  pixel  is  A^,  the  number  of  resolvable  bands 
L  is  limited  by  L  =  a  A2^Al  •  The  spectral  resolution  is  limited  by  Additionally,  it  is  assumed  that 
the  side  length  of  the  coded  aperture  square  pixel  Ac  satisfies  kAc  =  A^,  where  k  >  1  is  an  integer. 
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Figure  3.1:  Illustration  of  the  spectral  data  flow  in  CASSI.  The  qth  slice  of  the  data  cube  T  with  11 
spectral  components  is  coded  by  a  row  of  the  coded  aperture  t  and  dispersed  by  the  prism.  The  detector 
captures  the  intensity  y  by  integrating  the  coded  light. 


The  horizontal  and  vertical  spatial  resolutions  are  thus  limited  by  A^.  The  discretized  output  at  the 
detector  is  given  by 


L—l 


Gjk  ^  y  T  OJ jk 


(3.1) 


i= 0 


where  Gjk  is  the  intensity  at  the  j,  k  position  at  the  detector  whose  dimensions  are  N  x  (M  +  L  —  1). 
T  is  an  L  x  N  x  M  spectral  data  cube,  T  is  the  coded  aperture  and  uj  is  the  white  noise  of  the  sensing 
system.  The  output  Gjk  can  be  written  in  matrix  notation  as 


yt  =  $/  +  (jJ 


(3.2) 


where  yt  is  a  vector  representation  of  G  and  /  is  the  vector  representation  of  the  data  cube  T .  Notice  that 
/  is  represented  as  /  =  ^3 dP  where  p  are  the  sparse  coefficients  of  representation  on  the  basis  [?]. 
In  this  Chapter,  the  coded  aperture  is  considered  binary  and  the  dispersive  element  is  considered  linear. 
In  practice,  it  is  necessary  take  into  account  the  various  optical  artifacts  and  non-ideal  characteristic  of 
the  optical  system.  Observe  in  (3.1)  that  the  spatial  dimension  j  can  be  fixed  to  a  specific  value,  for 
example  j  =  q  represents  the  qth  row  measurement.  In  this  case  (3.1)  would  represent  the  model  of  a  slice 
of  CASSI.  Figure  3.1  shows  a  typical  slice  model  where  a  slice  F  of  the  data  cube  T  is  coded  by  the  coded 
aperture  t,  dispersed  by  the  prism,  and  the  resulting  light  is  integrated  at  the  detector  resulting  in  the 
measurement  y.  The  vectorization  of  F  is  f.  Given  the  separability  in  slices  of  (3.1),  a  slice  model  of  the 
instrument  is  developed  first.  The  results  are  thereafter  generalized  to  a  complete  2D  detector.  Observe 
that  the  coded  aperture  effects  are  embedded  in  the  matrix  $  in  (3.2),  then  a  more  suitable  matrix-based 
model  than  that  given  by  (3.2)  needs  to  be  developed  in  order  to  optimize  the  coded  aperture.  This  is 
developed  in  the  remainder  of  this  section.  To  begin,  let  0l  be  the  L  x  L  cyclic  permutation  matrix 


O/  = 


0  0  ...  0  1 

1  0  ...  0  0 

0  1  ...  0  0 

:  0  *  •  0  0 

0  0  ...  1  0 


(3.3) 


Next,  define  the  matrix  as  the  Lx  L  matrix  having  only  one  non-zero  element  in  the  coordinate 
as 


1  0  ...  0 

0  0  ...  0 


Jl  = 


_  0  0 

To  simplify  the  notation  define  the  matrix  Pl;/c  as 


0 

0 


(1,1) 


(3.4) 


pL;k  =GkLJL(Gl)k. 


(3.5) 
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Matrix  F 


Figure  3.2:  The  qth  slice  of  the  data  cube  T  is  represented  by  the  matrix  F.  Each  Fjk  element  is  pictorially 
represented  as  a  small  cube  where  the  gray  color  indicates  a  zero  value. 


Notice  that  multiple  applications  of  the  matrix  in  (3.3)  perform  the  same  cyclic  permutation  several 
times.  The  matrix  Pl;/c  in  (3.5)  will  be  used  to  select  an  specific  row  or  column  of  a  matrix.  Finally, 
define  the  reversing  order  matrix  C  as 


"  0 

0  .. 

..  1 

c  = 

0 

1 

'•  0 

(3.6) 

1 

0  .. 

..  0  _ 

3.2.1  Slice  Model  of  the  Spectral  Data  Cube 

The  qth  slice  of  T  is  created  as 

Fjk=Fqjk  j  =  0,  — 1,  k  =  0, . . .  ,  L  -  1.  (3.7) 

Notice  that  each  column  of  F  contains  spatial  information  at  a  given  wavelength  and  each  row  of  F 
consists  of  all  spectral  components  at  a  specific  spatial  position.  Figure  3.2  depicts  a  typical  matrix  F 
representing  a  slice  of  a  spectral  data  cube  and  the  respective  elements  present  in  CASSI.  The  kth  column 
of  F  is  denoted  as  f then  the  matrix  F  can  be  expressed  as  F  =  [f0,  fi, . . . ,  Observe  that  F  can 

be  decomposed  into  L  matrices  as 

L- 1 

F  =  E  Ffc  (3.8) 

k= 0 

k  N-l-k 

where  F k  =  [Om,  • . . ,  0 m  ,  ffc,  Om,  • . . ,  0 m\  and  0 m  is  an  M  long  column  vector  with  zero  elements. 
Observe  that  F k  can  be  written  more  succinctly  as 

F/c=FPl;/c  (3.9) 

where  P L-k  is  given  in  (3.5).  The  slice  F  can  then  be  rewritten  as 

L- 1 

F  =  ^  FPL;fe.  (3.10) 

k= 0 


3.2.2  Coded  Aperture  Effect  Model 

Consider  the  slice  F  in  (3.10)  impinging  on  a  row  of  the  2D  coded  aperture  T,  as  depicted  in  Fig.  3.1. 
This  row  of  the  coded  aperture,  represented  by  the  vector  t  =  [to, ... ,  Im- i]  ,  modulates  the  columns  of 
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F.  The  output  after  the  coded  aperture  can  be  represented  by  F  =  TF  where  T  =  diag  (to, . . . ,  £m-i)  is 
an  M  x  M  diagonal  matrix  containing  the  respective  row  of  the  coded  aperture.  As  mentioned  previously, 
the  coded  aperture  t  is  defined  as  the  Hadamard  product  of  two  vector  components  as  t  =  ?  o  w,  where 
?  =  [fo,  •  •  • ,  tm—i\  accounts  for  a  pseudorandom  component  and  w  =  [wq,  . . .  ,wm-i\  accounts  for  an 
adjustable  component  of  t.  T  can  thus  be  expressed  as  T  =  RW  where  R  =  diag(fo, . . . ,  tm-i)  and 
W  =  diag  (u)o, . . . ,  wm- i).  Figure  3.2  shows  that  in  CASSI  only  L  pixels  at  the  detector  are  affected  by  a 
point  source  in  the  spectral  data  cube.  More  specifically,  in  Fig.  3.2  the  component  of  the  coded  aperture 
(t)M-i  with  the  dashed  outline  affects  only  L  pixels  at  the  detector  y.  Hence,  both  the  adjustable 
component  w  and  the  random  component  ?  can  be  modeled  as  L-long  vectors  repeated  via  a  modulo 
operator  to  span  the  length  M  of  the  detector.  More  specifically  the  elements  of  w  and  ?  are  calculated 
as 


Wj  ^mod L(j)  3  1 

rj  =rmodL(j)  j  =  0,  •  •  • ,  M  -  1  (3.11) 

where  w  and  r  are  L-long  vectors.  Notice  that  the  coded  aperture  elements  can  be  expressed  as  tj  = 
fmod L(j)  for  j  =  0, . . . ,  M  —  1  where  the  vector  t  is  t  =  r  o  w.  These  vectors  can  thus  be  written  as 


W  =U M'  0  w 
r  =uM'  0  r 

t  =uMf  0  t  (3.12) 

where  0  is  the  Kronecker  product,  and  u m'  is  the  one- valued  vector  of  length  M'  =  From  (3.12) 

WandR  can  be  expressed  as  W  =  diag  (u m>  0  w)  and  R  =  diag  (u m'  0  r).  The  matrix  F,  representing 
the  coded  slice  F,  can  thus  be  expressed  as 


F  =  RWF. 


(3.13) 


3.2.3  Dispersive  Element  Effect  Model 


Given  the  matrix  representation  of  the  coded  slice  at  the  source  in  (3.13),  here  we  show  the  effect  of 
the  dispersive  element  over  this  matrix.  Observe  that  each  column  of  F  represents  a  given  spectral  band 
spatially  coded  by  the  coded  aperture.  The  effect  of  the  dispersive  element  can  then  be  modeled  as  a 
shift  of  each  column  of  F  in  k  units  where  k  expresses  the  kth  wavelength  (column).  Further,  the  matrix 
F  needs  to  be  augmented  by  L  —  1  rows  to  correctly  model  the  prism’s  effect.  Let  the  matrix  Y  model 
the  output  of  the  dispersive  element  when  the  input  is  the  matrix  F  =  [f0, . . . ,  ?l-i]  in  (3.13).  Y  is  then 
given  by 


'  fo 

Oi 

Y  = 

fi 

Ol-i 

Ol_i 

O 

1 

to 

- 1 

i 

|t4H 

where  0^  =  O.u^  is  an  i— long  zero- valued  column  vector.  Figure  3.2  illustrates  the  formation  of  Y.  Define 
V  as  F  =  M  +  L  —  1.  In  order  to  express  the  matrix  Y  compactly,  the  V  x  L  matrix  F&  is  defined  as 


Ovx/c 

f  k 

Ovx(L-fc-l) 

0(L— l)xl 

(3.15) 


where  Oyx/c  and  0yx(L_/c_1)  are  matrices  of  V  x  k  and  V  x  (L  —  k  —  1)  zero  elements  respectively.  Then 
the  matrix  Y  can  be  written  as 

L- 1 

Y=^(0y)feFfc  (3.16) 

k= 0 

1If  ^  is  not  an  integer  then  the  first  elements  of  w  and  ?  are  given  by  (3.12)  and  the  remaining  by,  Wj  =  Wj-L , 

and  fj  =  Tj-L  f°r  3  =  L  \_1^\  >  •  •  •  >  M  —  1. 
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where  the  V  xV  matrix  has  the  same  purpose  and  structure  of  0^  in  (3.3).  Defining  the  matrix  I  as 

1=  [ImxmOmx(l-i)]T  (3-17) 

where  I mxM  is  an  M  x  M  identity  matrix,  Fk  in  (3.15)  can  be  expressed  as 

Fk=IKWFk  (3.18) 

where  Fk  is  given  by  (3.9).  Replacing  (3.9)  and  (3.18)  in  (3.16)  the  matrix  representation  of  the  prism 
effect  on  F  is 

L- 1 

Y  =  ]T(ev,)felRWFPL;fe.  (3.19) 

k= 0 

Observe  that  the  jth  row  of  Y  contains  the  spectral  information  that  is  integrated  in  the  jth  pixel  detector. 


Table  3.1:  CASSI  Equations  Model  Summary. 


Dimensions 

Complementary  Equations 

Data  Cube 

Full  Data  Cube  T 

N  xM  x  L 

Fj  —  J~ qjki 

Slice  F  =  FPL;fc 

M  x  L 

F  =  [fo  •  ■  ■  fi  - 1  ]  ■ 

Coded  Aperture 

2D  Coded  Aperture  T 

N  x  M 

R  =  diag  (?),  W  =  diag  (w), 

Slice  Coded  Aperture  t 

L 

T  =  RW,  w  =  uM/  ®  w, 

Random  Component  r 

L 

f  =  uMf  O  r,  t  =  UM/  O  t, 

Designable  Component  w 

L 

t  =  r  o  w. 

Coded  Aperture  Output  F  =  RWF 

M  x  L 

Prism  Output 

V  =  M  +  L- 1 

Y  =  E^=n(©M)fclRWFPL;fc 

V  x  L 

VL;k=ekr3L{&l)k- 

Detector  Output  ,  two  versions: 

y  =  Epo1  Py;  E^01(©^)fclFPL;fcc(eD^«t 
y  =  Etc1  pv-,j  Etc l(&v)kiRFPL;kc(eiy+1w 

V 

V 

3.2.4  Detector  Measurements 

The  detector  integrates  all  coded  light  that  comes  from  the  dispersive  element.  The  vector  of  intensities 
at  the  detector  can  be  calculated  as 

y  =  YuL  (3.20) 

where  is  a  one- valued  L  long  vector.  Substituting  (3.19)  in  (3.20),  we  obtain 

L—l 

y  =  ^(0v)fcIRWFPi;feuL.  (3.21) 

k= 0 

Property  1  Given  a  slice  F  as  input  to  CASSI  impinging  onto  the  slice  aperture  code  t,  the  FPA 
measurement  at  the  detector  is 

v-i  L- 1 

y  =  51  pv-j  E (0 v)kiFPL.kc(eTLy+H .  (3.22) 

j= 0  k= 0 

Equivalently ,  the  output  is  given  by 

v-i  L- 1 

y  =  X]  pV;j  E  (6y )feIRFPL;fcC(eD^w,  (3.23) 

j= 0  k= 0 

where  P y.j,  0m,  I,  Pl;/c,  C,  and  0^  are  the  shifting  matrices  defined  in  (3.5),  (3.3),  (3.17),  and  (3.6). 

Observe  that  (3.22)  is  specially  suitable  for  the  design  of  the  coded  aperture  t  since  t  appears  at  the 
right  most  element  in  (3.22).  In  contrast,  the  most  right  element  in  (3.23)  is  w.  Hence,  Eq.  (3.23)  can 
be  used  to  design  the  component  w  of  the  coded  aperture  when  the  random  component  r  is  assumed  to 
remain  fixed.  Table  3.1  summarizes  the  principal  results  of  the  slice  CASSI  model. 


52 


(e) 

Figure  3.3:  (a)  A  random  coded  aperture  and  (b)  a  spectrally  selective  optimal  coded  aperture  for  a  12 
shot  CASSI  system.  The  corresponding  FPA  measurements  are  shown  in  (c)  and  (d).  The  zoomed  areas 
illustrate  the  wavelengths  present  at  each  pixel  measurement  where  the  spectral  selectivity  of  the  optimal 
codes  is  clearly  seen.  The  desired  spectral  profile  A  E  [461nm  —  471nm,  641nm  —  668nm]  is  illustrated 
in  (e). 

3.3  Optimal  Codes  for  Spectral  Band  Selectivity 

Spectral  selectivity  aims  at  designing  coded  aperture  sets  that  obtain  compressive  measurements  with 
only  information  from  a  given  set  of  bands  of  interest.  Furthermore,  the  optimization  aims  at  the  use  of 
the  smallest  number  of  measurements.  Figure  3.3  illustrates  a  sample  of  a  spectrally  selective  optimal 
coded  aperture  and  the  respective  FPA  measurement.  There,  it  is  possible  to  observe  that  the  FPA 
spectrally  selective  measurement  contains  only  wavelengths  from  a  specific  set  of  bands.  In  particular, 
observe  that  each  pixel  in  the  zoom  area  in  Fig  3.3(d)  contains  exclusively  spectral  components  in  the 
range  [461nm  —  471nm,  641nm  —  668nm].  In  contrast,  the  zoom  area  in  Fig.  3.3(c)  shows  that  a  CASSI 
measurement  using  random  codes  contains  spectral  components  of  the  whole  wavelength  range.  The 
coded  aperture  design  developed  hereafter  can  be  summarized  as  follows:  A)  Obtain  an  expression  for 
the  desired  spectrally  selective  compressive  measurement  signal  d  from  the  observation  spectral  image 
slice  F.  B)  Determine  a  set  of  K  optimal  weights  {wQj.  such  that  the  vector  g0  =  d  can  be 

constructed  from  the  respective  K  compressive  measurements  {yaj  •  C)  Generalize  the  results  in 
(3.3)  by  incorporating  the  random  components  of  the  coded  aperture  raj  such  that  the  set  {gj}^=Q 
containing  U  different  realizations  of  the  desired  signal  d  can  be  estimated.  D)  Optimize  the  realizations 
of  the  pseudorandom  components  raj  such  that  the  number  of  shots  K  in  (3.3)  is  further  reduced  to  a 
desired  number  K'  <  K. 


3.3.1  Desired  Compressive  Measurements 


The  reference  compressive  measurement  signal  d  is  created  by  zeroing  out  all  spectral  bands  at  the  input 
except  those  that  are  desired.  More  specifically,  define  the  L  long  vector  A  indexing  the  desired  spectral 
bands  as 


A  k  — 


if  the  kth  spectral  band  is  of  interest 
otherwise. 


(3.24) 


Furthermore,  let  C  =  ||A||o  be  the  number  of  desired  bands,  define  the  index  of  the  jth  nonzero  component 
of  A  as  ej,  and  let  A  =  diag  (Aq,  . . . ,  Xl-i)  be  the  diagonal  matrix  containing  the  desired  band  indices 
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in  (3.24).  Given  the  slice  F,  the  matrix  FA  contains  only  the  spectral  information  from  the  desired 
bands.  Fixing  the  random  component  r  and  making  w  =  u^,  y  in  (3.23)  results  in  the  output  containing 
the  information  from  the  spectral  data  cube  F  coded  by  the  random  component  r.  In  addition,  if  F  is 
replaced  by  FA  then  the  output  is  a  compressive  spectral  measurement  containing  only  information  of 
the  spectral  bands  of  interest  coded  by  r.  More  specifically,  using  (3.23)  the  desired  output  is 

V-l  L- 1 

P vj  (©v)felRFPL;fcA.  (3.25) 

j=0  k= 0 

Notice  that  d  in  (3.25)  is  created  by  sensing  the  desired  data  cube  FA  and  by  bypassing  the  effect  of  the 
coded  aperture  w.  The  next  sub-section  aims  at  obtaining  the  same  values  of  the  vector  d  by  sensing 
the  complete  data  cube  F  using  the  CASSI  system  and  by  varying  the  term  w  in  (3.23). 


3.3.2  Coded  Aperture  Design  for  w 


Given  the  expressions  for  y  and  d  in  (3.23)  and  (3.25),  respectively,  the  objective  in  the  design  of  w  in 
(3.23)  is  to  make  all  the  elements  of  y  and  d  equal,  where  the  random  components  of  the  coded  aperture 
r  are  not  yet  considered,  such  that  all  the  elements  of  r  are  set  to  one.  In  order  to  estimate  the  weights 
w,  the  output  in  (3.23)  is  equated  to  that  of  the  desired  response  in  (3.25)  resulting  in 

A  =  C(0lr+1w  jm  (3.26) 

Solving  (3.26)  for  w  results  in 

w  =  (0Ly+1C-1A  j  =  0,...,V-l  (3.27) 


where  it  can  be  observed  that  a  single  vector  w  cannot  satisfy  (3.27)  for  all  j.  Hence,  V  —  1  vectors  must 
be  used.  Let  these  vectors  be  Wj  for  j  =  0, . . . ,  V  —  1.  Let  yj  be  the  CASSI  output  in  (3.23)  when  the 
coded  aperture  w j  is  used,  then  the  element  by  element  equality  condition  (d )j  =  (y )j  can  be  written  as 


(d)i  =  (y;)i  3=  o . V-l.  (3.28) 

Since  (0j)J  =  (0^)J+L,  then  only  L  values  of  j  are  needed  to  be  taken  into  account  in  (3.27),  leading  to 


wj  =  (0L)i+1  C-XA 

W  J+mL  =  W  j 


j  =  0,...,L-l 
,  V 

rn  =  0, 1 . . . 


(3.29) 


Consequently,  (3.28)  is  modified  to 


(d)j  —  (y  modi,  (j) )  J  j  —  O’  •  •  •  5  ^  1-  (3.30) 

The  condition  in  (3.30)  can  then  be  satisfied  using  L  CASSI  shots  whose  outputs  are  given  by  yo, . . . ,  yz,-i- 
Observe  that  some  of  coded  apertures  in  (3.29)  can  be  expressed  as  a  linear  combination  of  the  others. 
Further,  since  w j  and  yj  are  linearly  related,  some  of  the  outputs  yj  can  also  be  expressed  as  the  linear 
combination  of  others.  It  is  thus  possible  to  calculate  d  in  (3.30)  using  K  <  L  shots.  To  estimate 
the  K  linearly  independent  weight  vectors,  the  coded  apertures  in  (3.29)  are  arranged  into  the  matrix 
-  [wo, . . . ,  wl-i]  or  equivalently 

Mw  =  [GlC-'X,  (0l)2C-xA,  . . . ,  (0i)LC-1A]  .  (3.31) 

Then,  the  minimum  number  of  shots  K  for  a  given  set  of  bands  of  interest  is  the  number  of  independent 
columns  of  Mw.  Hence, 

K  =  rank(Mw)  (3.32) 

where  rank(Mw)  is  the  rank  of  the  matrix  Mw.  The  K  linear  independent  columns  of  Mw  are  the  linear 
independent  weight  vectors  wak  for  k  =  0, . . . ,  K  —  1.  The  remaining  L  —  K  columns  of  Mw  can  be 
estimated  using  the  ensemble  of  vectors  W  =  [wao, . . . ,  wttK_1] .  The  jth  vector  w j  in  (3.29)  can  be 
expressed  as 

K- 1 

Wj  =  ^2  ckjV?ak  for  j  =  0, . . . ,  L  -  1  (3.33) 

k= 0 

where  the  coefficients  Ckj  are  the  elements  of  C  —  (WTW)  1  WTMw. 
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3.3.3  Spectrally  Selective  Image  Measurements 

Given  the  optimal  coded  aperture  set  W  =  [wao, . . .  and  the  corresponding  CASSI  measure¬ 

ments  of  a  scene  yao, . . . ,  yaK_1?  the  goal  is  to  re-order  these  measurements  and  select  the  entries  which 
contain  only  the  bands  of  interest.  It  turns  out  that  the  redundant  representation  of  the  weights  shown 
in  (3.33)  provides  a  simpler  structure  in  which  to  identify  the  entries  of  interest.  To  this  end,  the  set  of 
measurements  {yaj  is  expanded  into  the  set 


K- 1 

y j  =  T  Ckjyak  j  =  o, . . . ,  L  -  1  (3.34) 

k= 0 

where  (3.34)  uses  the  fact  that  there  exists  a  linear  relation  between  w j  and  yj.  Note  that  in  obtaining 

(3.34) ,  we  are  not  taking  additional  measurements,  we  are  only  combining  entries  of  the  existing  set  of 
measurements  {ya.}.  Once  the  expanded  measurement  vectors  {y j}  are  available,  the  elements  within 
this  set  that  contain  only  the  spectral  component  of  interest  must  be  identified.  Let  {g \j}f=$  be  the  K 
possible  measurement  vectors  that  satisfy  the  above  condition  obtained  by  selecting  elements  in  (3.34). 
Property  2  describes  the  structure  of  the  first  vector  go. 

Property  2  Let  Z  =  [yo, . . . ,  Yl-i]  be  the  ensemble  matrix  containing  the  expanded  measurements  vec¬ 
tors  in  (3.34).  The  vector  go  G  obtained  from  Z  as  (go)*  =  (Z )*modL(d  or  equivalently  (go)*  = 
(ymodL(i))i  for  i  —  0,  •  •  • ,  V  —  1  contains  measurements  with  spectral  component  of  interest  only ,  regard¬ 
less  of  the  value  taken  by  the  pseudorandom  component  r. 

Proof:  1  When  the  pseudorandom  component  is  taken  into  account  Eq.  (3.27)  is  rewritten  as 

A  o  r  =  r  o  C(©J)J+1w  j  #;  0, . . . ,  V  —  1.  (3.35) 

The  optimal  set  of  weights  when  the  pseudorandom  component  is  taken  into  account  is  given  by  solving 

(3.35) .  Notice  that  if  the  set  of  vectors  {w^Iq-1  solve  (3.27)  then  they  also  solve  (3.35).  Hence,  the 
optimal  set  of  weights  calculated  without  taking  into  account  the  random  components  can  be  used  to  solve 

(3.35).  Consequently,  equations  (3.29)  and  (3.30)  which  were  derived  from  (3.27)  are  valid  regardless  of 
the  values  taken  by  the  pseudorandom  components.  From  (3.30),  note  that  (d)o  =  (yo)o>  (d)i  =  (yi)i; 
and  so  on,  such  that  it  is  possible  to  construct  a  selective  measurement  using  the  set  {y*}.  Let  go  be  the 
selective  measurement  whose  elements  are  given  by 

go  =  [2/00?  2/11?  •  •  •  j2/(Lt4)(L-l)?2/(0)L?2/l(L+l)j 

•  •  •  5  ymodL(V-l)(V-l)]  (3.36) 

where  the  element  yji  is  yji  =  (y j)*.  Given  the  matrix  Z,  Eq.  (3.36)  can  be  seen  as  (go)*  =  (Z )*modL(i) 
for  i  =  0, . . . ,  V  —  1. 

Note  that  the  elements  in  go  are  the  diagonal  entries  in  Z  as  illustrated  in  Fig.  3.4(a).  It  is  important 
to  point  out  that  many  other  elements  in  Z  contain  spectral  components  of  interest,  in  addition  to  other 
spectral  components  that  are  not  of  interest.  Note  however  that  these  undesirable  spectral  components 
may  disappear,  in  some  cases,  when  the  pseudorandom  component  of  the  coded  aperture  is  taken  into 
account.  In  fact,  the  pseudorandom  structure  of  r  can  be  designed  to  improve  the  measurement  process  as 
described  in  Section  3.3.4.  Figure  3.4  illustrates  this  observation  where  Fig.  3.4(a)  depicts  the  elements  in 
Z  which  contain  spectral  components  of  interest,  regardless  of  the  value  taken  by  r.  Figure  3.4(b)  depicts 
the  elements  in  Z  which  contain  spectral  components  of  interest  only,  when  r  is  taken  into  account  and 
where  r  is  a  realization  of  a  Bernoulli  random  variable  with  parameter  p  =  0.6.  Note  that  the  elements  in 
Fig.  3.4(a)  form  a  subset  of  the  elements  in  Fig.  3.4(b).  Figure  3.4(c)  illustrates  the  elements  of  interest 
when  r  is  a  realization  of  a  Bernoulli  random  variable  with  parameter  p  =  0.3.  Note  that  as  p  decreases, 
the  number  of  elements  in  Z  containing  the  bands  of  interest  increases.  It  will  be  shown  shortly  that  an 
optimal  p  is  obtained  in  the  optimization  of  the  codes.  The  entries  in  Z  depicted  in  Fig.  3.4(a)  are  used 
to  form  the  vector  go-  Next,  a  procedure  to  select  the  entries  depicted  in  Fig.  3.4(b)-(c)  which  were  not 
already  used  in  go,  is  identified.  Define  the  set  of  vectors  Xu  as 

A  it  =  C(0j)*+1(w^  o  Yf)  (3.37) 
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(a)  (b)  (c) 

Figure  3.4:  Three  different  versions  of  the  V  x  L  matrix  Z  as  the  value  of  p  decreases.  The  elements 
available  to  construct  the  set  of  optimal  measurements  do, . . . ,  du  are  shown  in  white  squares  representing 
elements  containing  only  the  desired  bands,  (a)  The  pseudorandom  component  is  not  considered,  pm  1. 
The  pseudorandom  component  is  a  realization  of  a  Bernoulli  variable  with  parameter  (b)  p  =  0.6  and  (c) 
p  =  0.3.  In  this  example:  L  =  8,  V  =  280. 


Table  3.2:  Spectral  Selectivity  Equations  Summary. 


Size 

Complementary  Equations 

Vector  of  desired  bands,  A 

L 

Xek  =  1,  0  otherwise.  A  o  =  £ 

Optimal  weight  apertures 

w,  =  (0Ly+1  c-xA 

Matrix  of  weights 

L 

K  Linearly  independent  columns 
of  1VU :  waj 

Mw  =  [w0, . . .  ,wl_i] 

Lx  L 

K  =  rank(Mw) 

Optimized  measurements 

Si  =  Ef=0  diag(njfc)ffc 

V 

{tj  =  UM'  ®  (rj  0  W.Ojto1 

r Z  ,\L  —  1  .  r rek\U  —  l,C  —  l 

l  J  j=0,fc=0 

n  jk  =  (©v)£fcI(t'fe) 

Ensemble  of  measurements 

4  =  (eoe‘  lfefc 

[g^>--->Sw-l]T  =Hb 

V  x  U 

[fo,  -  -  - ,  fc-l]  =^b 

H=  [Hi,..., Hu- i]T 

.  ..diag(n  =  'Hj  b 

for  i  =  0, . . . ,  V  —  1  and  faO,...,L-l,  which  indicates  the  spectral  bands  present  in  the  (i,  £)th  entry 
of  the  matrix  Z.  Using  these  vectors,  the  set  of  spectrally  selective  measurements  {g j}f=\  is  obtained 
by  forming  the  vectors  g j  as 


(g jh 


Zu  pe{0,...,I-l}:AtfVA-A  =  0 
0  otherwise 


(3.38) 


where  0  is  a  L-long  zero  valued  vector  and  it  is  assumed  that  the  elements  of  Z  used  to  construct  g j1 
are  not  used  to  construct  g j2  for  ji  <  j2.  Notice  that  the  number  of  linearly  independent  vectors  in  the 
set  is  less  or  equal  than  K.  Define  the  selectable  parameter  U  as  the  number  of  the  spectrally 

selective  measurements.  Then,  the  set  of  selective  measurements  is  limited  to  {gjl^Tg1  where  U  <  K.  The 
set  of  spectrally  selective  measurements  {gj}^=Q  in  (3.36)  and  (3.38)  can  be  interpreted  as  the  output  of 
a  modified  CAS  SI  system  similar  to  that  in  (3.21),  with  the  difference  that  each  spectral  image  slice  at 
the  kth  wavelength,  is  coded  by  a  different  coded  aperture  tk.  This  process  is  different  than  that  of  the 
traditional  CAS  SI  in  (3.21)  where  each  is  coded  by  the  same  coded  aperture  t.  Property  3  describes 
the  above  characteristic  in  the  modified  CASSI. 

Property  3  Given  the  set  of  spectral  bands  {I/d^To ,  the  matrices  Qy  and  I  in  (3.3)  and  (3.17),  the  set 
of  weights  {wjljCg1  in  (3.33),  the  set  of  vectors  representing  realizations  of  a  Bernoulli  random 

variable ,  and  given  the  set  of  coded  apertures  {t j  =  um'  ®  (rj  owj)}^1,  the  set  of  selective  measurements 
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{gj}y= o  can  be  expressed  as 


L- 1 


g  j  =  dias(*i)ffc 

k= 0 


where  (tj)i  =  (t f)i  for  some  I  G  {0, . . . ,  L  —  1}. 


(3.39) 


Proof:  2  Consider  go  first.  Using  (3.21)  and  observing  that  FPl-^ul  =  f k,  a nd  RW  =  diag(t*),  the 
measurements  yj  can  be  written  as 


L- 1 

yf  =  ^(0v)feI  diag(*i)ffc  (3.40) 

k= 0 

for  j  =  0, ...  ,L  —  1.  Observe  that  (3.40)  is  t/ie  traditional  CASSI  model  where  each  spectral  band  zs 
coded  by  the  same  code  t j.  Using  (3.40),  the  equality  (go)*  =  (ymodL(i))i  * n  Property  2  can  be  written  as 


(go); 


(L- 1 


^(6y)feI  diag(tmodL(i))ffe 


V/c=0 


L—l 


^  ((0vOfcI  diag(tmodL(;))4)i 

k= 0 


L—l 

Y  (  diag(tmodL(;))ffc)i-fe 

fc=0 


L—l 

(^modL(o) i—k  ^k\-k  ' 

k= 0 


(3.41) 


Defining  the  L-long  vectors  t§  as  (t§)*  =  (tmodL(o)i-/e  for  &  =  0, . . . ,  L  —  1  and  i  =  0, . . . ,  V  —  1,  then 
(3.41)  can  6e  written  as  (go);  =  Efe=o  (*o);  (f*0;-/c  =  Efc=o  ((0^)feIf4  (to)*-  ^ence, 


L-l 

go  =  y^(0v)fcI  diag  (t§)  f/c.  (3.42) 

/c=0 

Note  that  each  spectral  band  in  (3.42)  zs  coded  by  a  different  equivalent  coded  aperture  t§.  Similar 
expressions  to  that  in  (3.42)  can  6e  found  for  gj  for  j  >  0  szzc/z  that  gj  can  be  written  in  terms  of  tk- .  In 
these  cases,  the  vectors  ik-  depend  on  the  set  of  random  components  {rj}^~Q  and  thus  a  closed  expression 
as  (3.42)  is  not  available,  however,  they  can  be  computed  using  (3.38). 


The  expression  in  (3.39)  can  be  further  simplified  observing  that  the  vector  gj  only  contains  information 
from  the  C  desired  bands  defined  in  (3.24)  and  indexed  by  the  vector  6j.  Thus,  the  vectors  gj  can  be 
rewritten  as 

c- 1 

Si  =  5T0vd£,cI  diag {tf%k  (3.43) 

k= 0 


for  j  =  0, . . . ,  K  —  1.  Define  the  vectors  =  (0y)efeI(t^fc)  and 
be  succinctly  represented  as 


c- 1 


g j  =  Y  diaS (njk)h 

k= 0 


A 


(0y)e/Tfefc,  then  gj  in  (3.43)  can 


(3.44) 


for  j  =  0, . . .  ,U  —  1.  Notice  that  and  in  (3.44)  are  circular  shifted  versions  of  the  vectors  te-k  and 


f£fc  respectively.  The  sparse  representation  of  the  ensemble  of  vectors  f  : 
the  vectors  gj  can  be  expressed  as 

g j  =  Hjb 


,  I£-1 


is  f  =  \I/b.  Then, 
(3.45) 
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where  /Hj=J2k=o(l1’  >  °]0l  )  ®  diag(nJ;c)'J>.  The  ensemble  of  the  vectors  g,  can  be  expressed  as 


n 


s'  ’  \ 


go 

H0 

gw-l 

_  Hu- 1  . 

b. 


(3.46) 


3.3.4  Design  of  the  Pseudo-Random  Components 

The  minimum  number  of  shots  K  given  by  (3.32)  assumes  that  the  pseudorandom  component  Yj  of  the 
coded  aperture  tJ  =  w j  o  Yj  remains  fixed.  It  is  shown  here  that  the  number  of  shots  can  be  further 
reduced  if  the  pseudorandom  component  is  optimized.  The  optimization  is  accomplished  as  follows.  The 
Hadamard  multiplication  of  both  sides  of  (3.29)  by  Yj  produces 

ri  °  wi  =  ri  °  ((©i)i+1C_1A)  (3.47) 

tj  =  R.,(0,.)-'flC  'A  (3.48) 

for  j  =  0, ...  yL  —  1,  where  R j  =  diagjrj}.  Using  (3.48),  the  matrix  M.w  in  (3.31)  can  be  generalized  to 
include  the  pseudorandom  components  as 

M t=  [to  ...  t^-i] 

=  [Ro0LC-1A,R1(eL)2C-1A,...,RL_1(0L)LC-1A]  . 

The  objective  of  the  optimization  of  the  set  of  pseudorandom  vectors  {1*0, . . . ,  is  to  minimize  the 

rank  of  until  a  predetermined  number  of  shots  Kr  <  K  is  obtained.  Some  constraints  are  necessary  in 
order  to  limit  the  search  of  the  vectors  Yj.  In  our  approach,  we  seek  the  vectors  Yj  such  that  the  CAS  SI 
measurements  better  satisfy  the  Restricted  Isometry  Property  (RIP)  condition.  The  RIP  is  an  important 
condition  for  compressive  measurements  which  guarantees  the  robust  recovery  of  the  sparse  data  cube 
via  1 1  minimization  [?,  ?,  ?]. 


Figure  3.5:  Optimization  of  the  coded  aperture  process.  Given  the  vector  A,  the  optimization  reduces 
the  rank  of  the  matrix  containing  the  set  {w j  o  Yj}^~^  where  the  varying  terms  are  the  vectors  Yj. 
The  optimization  is  constrained  to  satisfy  the  condition  given  in  (3.50). 


Property  4  Given  a  set  of  spectral  bands  of  interest  f  = 


f0, . . . ,  ic-i 


their  sparse  representation 


f  =  tfb,  assume  that  |b|  =  S ,  given  the  matrix  TL  defined  in  (3.46),  the  vectors  Ujk  defined  in  (3.44), 
and  the  constants  Sc  G  (0, 1)  and  e  G  [0.5, 1],  then  the  probability 


P((l-4)  II  b  ||2<||  Hb  \\l<  (1  +  5C)  ||  b  HI)  =l-e 


is  higher  when  the  constraints  ||  njk  llo  >M,  and  V  V  Y.%0  nTfcn^  =  0  for  M  =  0, . . . ,  C- 1 

are  satisfied.  Thus ,  if  a  CASSI  model  matrix  TL  satisfied  the  above  condition  with  probability  1  —  e\  and 
constant  Sc  then  that  matrix  can  be  optimized  to  satisfy  the  mentioned  constraints  such  that  the  above 
condition  is  satisfied  for  the  same  constant  Sc  with  probability  1  —  62  where  e\  >  e 2. 
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The  constraint  ||  X^=oni&llo  >  M  accounts  for  the  necessary  condition  to  reconstruct  each  spectral 
band  .  The  constraint  ££fc,  nJknjt  =  0  accounts  for  the  joint  reconstruction  of  the  complete 

set  of  desired  bands  in  the  data  cube  f .  These  two  conditions  are  used  as  constraints  in  the  optimization 
of  the  pseudorandom  terms  of  the  coded  apertures. 

Note  that  the  ith  element  of  the  vector  n  =  Y^=o  *22 f=o  nj 't  represents  the  approximate  number  of 
spectral  bands  sensed  at  the  ith  location  at  the  detector.  This  approximation  comes  from  the  correlation 
that  exists  among  the  vectors  n ju  and  Uj2£  for  all  t  and  j\  ^  Hence,  the  third  constraint  in  the 
optimization  is  ru  ^  0  for  i  =  0,  ...,H  —  1.  The  fourth  constraint  in  the  optimization  is  the  mean 
number  of  spectral  components  present  in  each  pixel  measurement  at  the  detector.  The  parameter 
/i*,  referred  as  the  optimal  compressive  ratio  is  a  selectable  term  that  indicates  the  average  number 
of  spectral  components  in  each  pixel  measurement,  and  thus  /jl*  >  1.  The  maximum  mean  number  of 
spectral  component  by  pixel  detector  can  be  estimated  as 

Mm  =  (3.49) 

where  ML  is  the  number  of  total  variables  and  V  is  the  number  of  effective  pixel  sensors  in  a  line  of  the 
FPA.  The  optimal  compressive  ratio  parameter  //*  range  is  thus  [1  In  general,  the  mean  number  of 
spectral  components  by  pixel  can  be  estimated  as:  y^Uyii,  where  uy  is  a  one- valued  V  long  vector. 


Table  3.3:  Iterative  Stochastic  Algorithm  to  solve  (3.50). 


Inputs 

A,  K/ ,  U ,  Pr(.)  is  a  random  permutation,  rk(.)  =  rank(.). 

Initialize 

w i  =  (0Lb+1C-1A,  M  =  [w0, . . .  M  =  M, 

B  =  {Bkj  =  Bernoulli(0.5)  :  J2k-o  Bkj  =  1  Vj, 

Efca  II  Efjo1  Bkj  ||o=  Kf}. 

Loop-1 

While,  rk(Ad)  ^  K' ,  A4  =  MoB,  return  to  Main. 

Loop-2A 

S  =  {Sj  e  {0, . . . ,  L  -  1}  :  Zt=o  M* ji  =  °}>  z=Pr(s). 

For  j  =  1, . . . ,  z  ;  for  (r,  m)  =  1,...,L;  r/m 

M'z.r  =  1,  M'z.m  =  1, 

If  rk(AT)  =K',M  =  M',  Co  =  0  break; 

Else  M'  =  M. 

return  to  Main. 

Loop-2B 

s  =  {Sj  e  {0, . . . ,  L  —  1}  :  Et=o  =  °}>  z=Pr(s). 

For  j  =  0, . . . ,  |z  —  1 1 ;  for  (r,  m)  =  0, . . . ,  L  —  1;  r  ^  m 

M'rZj  =  1,  M'mzj  =  1, 

If  rk(Al/)  ==  K' ,  M  =  A4/ .  ^i=0  break; 

Else  M'  =  M. 

return  to  Main. 

Loop-2C 

For  0, . . . ,  L  —  1;  x,  y  =  UniformDiscrete(0,  L  —  1) 

M’xy  =  1, 

If  rk(Ai')  =  K' ,  M.  =  M! ,  £2  =  0,  break; 

Else  M'  =  M. 

return  to  Main. 

Main 

While,  Si<M  V52> 0V  ||n||0<  V  V  /  >  fA  Vrk(WTW)  +  K' 

£  =  U3,  Loop-2  A,  Loop-2B,  Loop-2C; 

Update  nik,  n,  /  =  \/i*  -  ^u^n|, 

Wkm  =  MkoCrn  k  =  0, . . . ,  L  —  1,  m  =  0, .  . . ,  K'  —  1. 

If  £  || q=  3  then  Loop-1. 

Output 

Optimal  M*  =  [tg  . .  ,  W*  =[t*0, . . . ,  . 

c  =  (W*TW*)_1  W*M*. 

Given  the  vector  A,  the  respective  set  of  weights  {w j}f=g,  the  desired  number  of  shots  A',  the 
desired  number  of  optimal  measurements  U,  the  optimal  compressive  ratio  /i*,  the  optimization  of  the 


59 


pseudorandom  component  of  the  coded  apertures  r0, , . . ,  r^_i  can  be  written  as 


min  |/i* 

{r0v)rL-i} 


vu 


Subject  to 

rank([w0or0  ...  wL_i  o  rL_i])  =  K' , 
u-i 

\\J2n^\\o>M,  || n||0  =  V 

3=0 

k^t  j 


(3.50) 


where  the  vectors  representing  shifted  versions  of  the  coded  apertures  are  given  in  (3.44),  and 
n  =  Y^=o  represents  the  mean  number  of  spectral  bands  sensed  at  the  jth  location  at  the 

detector.  A  good  range  of  / 1 *  has  been  empirically  obtained  as  [1.4,  2.0].  In  addition  to  the  constraints  in 

(3.50) ,  the  vectors  r0, . . . ,  i\l_i  are  realizations  of  a  random  variable  perturbed  such  that  they  satisfy  the 
constraints  in  (3.50).  Figure  3.5  shows  a  flowchart  of  the  optimization  of  the  coded  apertures.  Solving 
equation  (3.50)  is  an  NP  hard  optimization  problem.  A  stochastic  based  algorithm  is  proposed  to  solve 

(3.50)  [?]. 


The  objective  of  the  algorithm  in  Table  3.3  is  to  find  the  matrix  Af  *  =  [tg  . . . ,  t x'-i]  which  satisfies 
rank(A4*)  =  K7,  and  the  constraints  given  in  (3.50).  Further,  the  optimal  matrix  Ai*  can  be  factor¬ 


ized  as  a  function  of  its  K'  linear  independent  columns  W*  = 


j.  *  j.  * 

bn  5  *  *  *  5  Lo 


as  M*  =  W*C  where 


e  =  (w*Tw*)  1  wtTM\  The  matrix  W*  contains  an  optimal  row  of  coded  apertures  for  spectral 
selectivity. 


The  algorithm  uses  the  LxL  matrix  B  whose  elements  are  random  realizations  of  a  Bernoulli  variable, 
further,  K'  rows  of  B  have  at  least  a  nonzero  element  or  the  same  rank(B)  <  K7.  The  Loop-1  uses  several 
realizations  of  the  matrix  B  to  initialize  the  matrix  Ai.  The  initial  matrix  Ai  has  both  rank(Af)  =  K7 
and  a  very  low  Frobenius  norm.  In  each  iteration,  the  Main  loop  verifies  that  the  constraints  given 
in  (3.50)  are  satisfied.  Then,  the  Main  loop  applies  the  loops  Loop-2A,  Loop-2B,  and  Loop-2C  to 
increase  the  Frobenius  norm  of  the  initial  matrix  Ai.  The  loop  Loop-2  A  increases  the  Frobenius  norm 
of  Ai  by  making  two  randomly  chosen  elements  in  a  column  of  Ai  equal  to  one.  The  loop  Loop-2B 
increases  the  Frobenius  norm  of  Ai  by  making  two  randomly  chosen  elements  in  a  row  of  Ai  equal  to  one. 
The  loop  Loop-2C  increases  the  Frobenius  norm  by  placing  a  one  element  in  an  arbitrary  coordinate  of 
Ai.  The  loops  Loop-2  A,  Loop-2B,  and  Loop-2C  do  not  change  the  rank  of  Ai.  Further,  the  Main 
loop  verifies  that  rank(WTW)  =  K7.  The  Frobenius  norm  of  Ai  is  increased  until  that  /  is  less  than 
some  user  defined  constant  /a-  If  the  loops  Loop-2  A,  Loop-2B,  and  Loop-2C  fail  simultaneously  to 
increase  the  Frobenius  norm  of  Ai,  then  the  algorithm  restarts  at  Loop-1.  Notice  that  the  algorithm 
in  Table  3.3  needs  to  be  run  N  times  to  construct  the  complete  2-dimensional  (2D)  coded  apertures.  A 
flowchart  of  the  procedure  to  reconstruct  the  desired  set  of  bands  is  shown  in  Fig.  3.6.  The  details  of  the 
construction  of  the  2D  coded  apertures  are  explained  in  the  next  section. 


3D  scene  fuk 


Figure  3.6:  Given  K'  2D  optimized  coded  apertures,  K'  CASSI  measurements  Yaj  are  made.  These 
measurements  are  reordering  to  construct  the  spectrally  selective  measurements  Gj.  The  GPSR  algorithm 
reconstructs  only  the  desired  bands.  N  unidimensional  optimal  codes  {rai  wttj.  }jLo  1  are  use(^  1°  construct 
the  2D  codes 
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Table  3.4:  Comparison  Between  CASSI  and  Optimized  CASSI 


CASSI 

Optimized  CASSI 

Number  of  Shots 

K 

K’  <K 

Data  Cube  Size 

N  x  M  x  L 

N  x  M  x  C 

Coded  Aperture 
for  each  mth  shot 

Constant,  Tm 

Variable  in  each  ith 
band  Tmi 

Output,  Gnk  = 

y  Fi,n,k+iTn,k+i 
i= o 

c-i 

\  "  ir>  rpi 

/  v  J  n,fe+z 

i= 0 

Measurements  Size 

N  x  (M  +  L-l)  x  K 

N  x  (M  +  C  -  1)U 

3.4  Optimal  2-Dimensional  Coded  Apertures 


Section  III  provided  the  coded  aperture  design  for  a  single  slice  at  the  detector.  Here,  we  extend  these 
concepts  to  the  design  of  2D  coded  apertures  which  are  built  by  an  ensemble  of  ID  coded  apertures.  The 
nth  execution  of  the  algorithm  in  Table  3.3  provides  the  optimal  vectors  , . . . ,  t^, .  The  N  x  K' 

set  of  vectors  {t^?}  are  used  to  create  the  K'  set  of  2D  coded  apertures  Tak  which  can  be  calculated  by 


(TaJni  =  (<&>) 


mod  L(j) 


(3.51) 


for  n  0, . . . ,  N  —  1,  j  =  0, . . . ,  V  —  1,  and  k  =  0, . . . ,  K'  —  1.  Let  Gak  be  the  2D  CASSI  output  when  the 
respective  coded  aperture  Tak  in  (3.51)  is  used.  Then,  Yak  is  the  2D  version  of  yak  in  (3.34).  Similarly, 
the  matrices  Y&  representing  the  2D  version  of  y*,  in  (3.34)  can  be  estimated  as 


K'-l 

Y k  =  ^2  Bn»fcYam  k  =  0,  .  .  .  ,  K  -  1 

m= 0 


(3.52) 


where  B 


mk  — 


diag  (c^k, . . . ,  C^k  ^ ,  and  the  matrices  C ^  are  provided  by  the  optimization  algorithm. 


Let  the  matrices  Y&  in  (3.52)  be  expressed  in  terms  of  their  columns  as  Y^  =  [y§  . . .  Yv-i]  ?  then  the 
2D  version  Go  of  the  desired  compressive  measurement  go  in  (3.36)  containing  only  the  bands  of  interest 
is  given  by 


G0  = 


y°, 


vL_1  v° 

’  y  l— i  ’  y  Li 


mod^CV- 

lYy-l 


(3.53) 


The  remaining  optimized  matrices  Gi, . . . ,  G^_i  are  estimated  using  (3.38)  and  selecting  the  elements  of 
the  matrices  Y&  in  (3.52)  containing  exclusively  the  desired  bands.  Observe  that  the  spectrally  selective 
measurement  matrices  Gm  in  (3.53)  can  be  seen  as  the  output  of  a  multi-frame  CASSI  system  sensing 
an  N  x  M  x  C  spectral  data  cube  given  by 


F'jk^Adk  i  =  0,...,£  (3.54) 

where  indexes  the  desired  bands.  Further,  the  matrices  Gm  can  be  seen  as  the  output  of  the  system 
whose  input-output  relation  is  given  by 

c- 1 

(Gm)ife  =  E  K^+iT^+i  +  ^,fe  (3-55) 

i= 0 


where  Trni  = 


fei(°)  fed1) 


Za(N-l) 


1  T 


and  tm^  is  the  nth  realization  of  the  vector  defined  in 

(3.43).  Table  (3.4)  summarizes  the  principal  differences  between  a  CASSI  system  using  multiple  shots 
and  the  CASSI  system  with  optimized  coded  apertures  defined  in  (3.55). 


3.5  Simulations 

To  test  the  methodology  developed  in  Sections  3.3,  and  3.4,  a  set  of  optimal  aperture  codes  for  two 
different  desired  spectral  band  sets  A  are  derived,  for  an  increasing  number  of  measurement  shots.  The 
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Figure  3.7:  4  spectral  channels  of  24  channels  in  the  data  cube  used  in  the  simulations  are  presented. 
The  complete  data  cube  extends  from  460nm  to  680nm,  and  it  has  24  spectral  channels  and  512  x  512 
pixels  of  spatial  resolution. 


experiments  use  a  test  data  cube  T  with  512  x  512  pixels  of  spatial  resolution  and  L  —  24  spectral  bands 
ranging  from  460nm  to  668nm.2  Figure  3.7  shows  a  portion  of  the  data  cube  used  in  the  simulations. 
The  compressive  sensing  reconstruction  is  realized  using  the  GPSR  algorithm  [?].  Sparsity  basis  selection 
for  hyperspectral  signals  has  been  widely  studied  in  [?].  It  is  shown  there  that  the  Kronecker  product 
of  ID  sparsity  basis  provides  the  highest  compression  ratio  for  hypersectral  images.  Here,  the  base 
representation  is  the  Kronecker  product  of  three  basis  T  =  Ti®4'2G)4'3,  where  the  combination 
Ti  (8)  T2  is  the  2D- Wavelet  Symmlet  8  basis  and  T3  is  the  Cosine  basis.  Other  Kronecker  products  such 
as  tridimensional  wavelets  can  also  be  used,  however,  they  requires  that  the  underlying  signal  has  dyadic 
dimensions.  The  basis  based  on  the  Cosine  Transform  has  the  advantage  of  its  low  computational  cost 
compared  with  the  tridimensional  wavelet.  Using  the  Kronecker  representation  the  data  cube  in  Fig. 
3.7  can  be  represented  with  only  4%  of  the  coefficients  and  preserving  the  99.67%  of  the  energy  of  the 
signal.  The  maximum  number  of  iterations  in  the  GPSR  algorithm  is  50  and  the  stopping  criterion  as 
recommended  in  [?]  is  the  objective  function. 


3.5.1  Rank  Minimization  Algorithm 

In  the  first  experiment,  the  rank  minimization  algorithm  is  run  for  L  =  24,  A  =  [461nm  —  471nm, 
641nm  —  668nm],  U  =  3,  K'  =  12,  /i*  =  1.6.  The  minimization  algorithm  aims  at  reducing  the  rank 
of  the  matrix  from  24  to  12.  Figure  3.8(b)  shows  the  error  /  =  |/z*  —  y^Uyii\  of  the  minimization 
algorithm  through  the  iterations.  There,  it  can  be  observed  that  the  algorithm  periodically  aims  at 
minimizing  the  cost  function,  however,  the  algorithm  restarts  when  the  constraints  in  (3.50)  are  not 
satisfied.  This  behavior  of  the  algorithm  is  originated  in  the  pseudorandom  searching  strategy. 

In  the  second  experiment,  the  optimization  algorithm  is  run  for  several  values  of  the  parameter  /i*, 
the  other  parameters  are  kept  equal  to  the  first  experiment.  For  each  /i*  value  the  optimization  algorithm 
is  run  512  times.  Using  the  spectral  data  cube  in  Fig.  3.7  and  the  512  realizations  of  the  optimization 
algorithm,  the  2D  spectrally  selective  measurements  {Go,Gi,G2}  are  constructed  using  the  procedure 
in  Section  3.4.  Using  the  spectrally  selective  measurements  and  the  GPSR  algorithm,  the  desired  bands 
are  reconstructed  and  their  mean  PSNR  is  calculated.  This  experiment  is  repeated  50  times  and  the 
PSNR  is  averaged.  Figure  3.8(a)  shows  the  averaged  PSNR  as  a  function  of  the  parameter  /i*.  From 

2Data  provided  by  Yuehao  Wu  and  Iftekhar  O.  Mirza.  Test  objects  provided  by  Andrew  Arce. 
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(a)  (b) 

Figure  3.8:  (a)  Performance  of  the  rank  minimization  algorithm  versus  the  /i*  parameter,  (b)  a  typical 
performance  profile  through  iterations  of  the  rank  minimization  algorithm.  The  iterations  where  the 
algorithm  restarting  are  indicated  with  squares. 


several  simulations,  it  was  found  that  regardless  the  value  of  L  or  C  the  optimal  value  of  the  parameter 
/i*  can  be  chosen  as  rule  of  thumb  between  1.4  to  2.0.  In  practice,  the  lower  the  value  of  /i*,  the  more 
difficult  it  will  be  to  achieve  convergence  of  the  optimization  algorithm. 


3.5.2  Spectral  Selectivity 

In  the  first  experiment  of  spectral  selectivity,  the  set  of  desired  bands  are  set  to  Ai  =  [461nm  —  479nm, 
641nm  —  668nm]  as  it  is  depicted  in  Fig.  3.9(g).  Figure  3.7  depicts  a  portion  of  the  24  x  512  x  512  data 
cube  used  in  the  simulations.  In  this  experiment,  the  initial  rank  of  the  matrix  is  24.  Using  the 
optimization  algorithm,  the  rank  of  the  matrix  is  minimized  to  9,  12,  and  18.  The  optimization  uses 
U  —  3  and  /i*  =  1.6.  The  experiment  is  repeated  45  times  and  the  mean  PSNR  is  estimated.  The  resulting 
spectral  data  cubes  are  shown  as  they  would  be  viewed  by  a  Stingray  F-033C  CCD  Color  Camera.  Figure 
3.9(a)  shows  the  desired  original  bands.  Figure  3.9(b)  illustrates  the  results  when  the  desired  spectral 
bands  are  reconstructed  using  12  Bernoulli  non-optimized  random  coded  apertures  measurements.  In 
this  case,  the  24  bands  are  reconstructed  and  the  desired  bands  are  manually  selected.  Figures  3.9(c)-(e) 
show  the  results  when  9,  12  and  18  optimized  coded  apertures  are  respectively  used  in  the  sensing  and 
reconstruction  process.  Notice  that  in  this  case  only  the  desired  bands  are  reconstructed.  Figure  3.9(f) 
shows  a  typical  realization  of  the  optimal  coded  aperture  mentioned  in  (3.55).  Figures  3.9(g),  3.9(h),  and 
3.9(i)  show  a  zooming  of  the  Figs.  3.9(a),  3.9(b),  and  3.9(d)  respectively.  Clearly  Fig.  3.9  depicts  that 
the  optimal  coded  apertures  conduce  to  a  higher  quality  reconstruction.  Additionally,  Figure  3.10  depicts 
the  differences  between  the  original  and  the  reconstructed  3rd  spectra  channel  (479nra).  Reconstructions 
artifacts  in  Figure  3.10  emerge  in  part  from  the  cosine  sparsifier  in  the  spectral  axis  and  in  part  from  the 
2D  wavelet  spatial  sparcifier. 


Table  3.5:  Mean  reconstruction  PSNR  in  dB. 


Number  of  shots 

9 

12 

18 

Ai 

A2 

Ai 

A2 

Ai 

A2 

Bernoulli,  min 

25.02 

26.25 

26.08 

27.13 

27.32 

28.12 

mean 

26.68 

27.47 

26.92 

27.77 

27.51 

28.30 

max 

27.19 

27.89 

27.38 

28.14 

27.58 

28.38 

Hadamard,  min 

23.51 

25.09 

24.66 

26.03 

25.71 

26.98 

mean 

25.50 

26.53 

26.15 

27.13 

26.99 

27.85 

max 

26.68 

27.37 

27.04 

27.83 

27.38 

28.19 

Optimized,  min 

30.79 

29.85 

29.44 

30.34 

31.75 

30.92 

mean 

30.91 

30.09 

31.02 

30.81 

31.79 

31.21 

max 

31.01 

30.27 

31.13 

30.96 

31.84 

31.32 

In  the  second  experiment,  9  desired  bands  are  selected  in  the  range  given  by  A2  =  [488nm  —  515nm, 
560nm  —  596nm].  In  this  case,  the  initial  rank  of  the  matrix  Mt  is  21.  The  optimization  uses  the  same 
values  ofW  and  / 1 *  used  in  the  first  experiment.  Again,  the  experiment  is  repeated  45  times  and  the  PSNR 
are  averaged.  Figure  3.11(a)  shows  the  desired  original  bands.  Figure  3.11(b)  illustrates  the  results  for 
12  Bernoulli  non-optimized  random  coded  apertures.  Figure  3.11(c)  shows  the  results  when  12  optimized 
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(j) 

Figure  3.9:  The  resulting  spectral  data  cubes  are  shown  as  they  would  be  viewed  by  a  Stingray  F-033C 
CCD  Color  Camera.  The  desired  bands  are  depicted  in  (g).  The  original  desired  bands  are  shown  in  (a). 
Reconstructed  images  by  using  :  (b)  12  shots  with  random  codes,  (c)  9,  (d)  12,  and  (e)  18  shots  with 
optimized  codes.  An  optimized  coded  aperture  used  to  reconstruct  (c)  is  shown  in  (f).  Zooming  of  (a) 
(original),  (b)  (random  codes),  and  (d)  (optimal  codes)  are  shown  in  (g),  (h),  and  (i)  respectively. 


coded  apertures  are  used  in  the  sensing  and  reconstruction  process.  Figure  3.11(d)  illustrates  a  typical 
realization  of  the  optimal  codes.  The  simulations  using  random  coded  apertures  were  repeated  using 
Hadamard  matrices,  however  the  PSNR  was  similar  to  that  of  the  Bernoulli  case.  Table  3.5  summarizes 
the  numerical  results  shown  in  Figs.  3.9  and  3.11.  To  indicate  the  consistency  of  the  algorithm  Table 
3.5  shows  the  minimum,  mean  and  maximum  averaged  PSNR  of  the  45  repetitions  of  the  experiment. 
The  optimized  results  indicate  improvements  up  to  4dB  in  the  PSNR  of  the  reconstructed  data  cubes 
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Figure  3.10:  Differences  between  the  original  and  the  reconstructed  3rd  spectra  channel  (479nra)  are 
presented  for  (a)  the  reconstructed  data  cube  in  Fig.  3.9(b)  (random  codes)  and  (b)  the  reconstructed 
data  cube  in  Fig.  3.9(d)  (optimized  codes).  A  zoomed  region  of  (a)  and  (b)  are  presented  in  (c)  and  (d) 
respectively. 


compared  with  the  results  of  using  random  coded  apertures. 


(e) 


Figure  3.11:  The  resulting  spectral  data  cubes  are  shown  as  they  would  be  viewed  by  a  Stingray  F-033C 
CCD  Color  Camera.  The  desired  bands  are  indicated  in  (e).  The  original  desired  bands  are  shown  in 
(a).  Reconstructed  images  by  using  :  (b)  12  shots  with  random  codes,  and  (c)  12  shots  with  optimized 
codes.  An  optimal  coded  aperture  is  illustrated  in  (d). 
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3.6  Conclusion 


A  matrix  representation  of  CASSI  has  been  developed.  The  matrix  model  permits  the  optimal  design 
of  weight  coded  apertures  that  can  be  used  to  build  a  set  the  spectrally  selective  compressive  CASSI 
measurements.  A  rank  minimization  algorithm  was  designed  to  estimate  the  pseudorandom  component 
of  the  coded  apertures.  The  optimization  of  the  coded  apertures  is  subject  to  a  specified  number  of  shots 
and  the  optimal  coded  apertures  are  designed  to  satisfy  the  RIP  condition  with  high  probability.  Given 
a  set  the  desired  bands,  the  PSNR  achieved  by  the  optimized  codes  is  up  to  4dB  higher  than  the  systems 
using  random  codes. 
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Chapter  4 


Higher-Order  Computational  Model 
for  Coded  Aperture  Spectral 
Imaging 

4.1  Introduction 

Typically,  the  discretized  output  at  the  CAS  SI  detector  gmn  is  modeled  as  the  sum  of  the  underlying 
spectral  voxel  slices  which  have  been  previously  modulated  by  a  coded  aperture  and  subsequently  spatially 
dispersed  by  a  prism.  More  specifically,  the  two  dimensional  CAS  SI  output  has  been  traditionally  modeled 
as 

L- 1 

Qmn  ^  ^  f  \m—k)jk'-^'(m—k)j  •>  (4*1) 

k= 0 

where  f\jk  is  the  discretization  of  the  underlying  spatio-spectral  power  source  density  and  is  the 
discretized  coded  aperture.  Notice  that  fijk  in  (4.1)  is  assumed  to  be  a  cubic  voxel  which  impinges  on  a 
single  pixel  detector  element  gmn.  The  analog  sensing  phenomena,  however,  is  such  that  when  a  single 
voxel  of  a  scene  is  dispersed  by  the  prism,  it  impinges  on  several  detector  elements  at  a  time.  This,  in 
turn,  causes  blurring  which  deteriorates  the  quality  of  image  reconstruction.  To  ameliorate  this  problem, 
instead  of  using  the  coded  aperture  T^-,  a  set  of  calibrated  coded  apertures  q  are  experimentally 

measured  and  used  in  the  reconstruction  process  to  take  into  account  the  non-ideal  optical  blur  and  non¬ 
linear  dispersion  [?].  Thus,  the  model  in  (4.1)  suffers  of  a  coarse  approximation  which  is  then  partially 
rectified  by  the  coded  aperture  calibration  process. 

There  are  two  drawbacks  in  this  approach.  First  the  calibration  process  is  often  inadequate  such  that  the 
discretized  voxels  f\jk  are  incorrectly  weighted  by  the  calibrated  codes  T^.  The  calibration  errors  orig¬ 
inate  principally  from  the  assumption  that  a  coded  cubic  voxel  impinges  on  the  detector  when  actually 
it  is  a  coded  oblique  voxel  which  impinges  on  it.  Secondly,  calibration  of  the  coded  apertures  is  difficult 
for  multiframe  CAS  SI  systems  where  a  sequence  of  coded  apertures  are  used  sequentially.  This  Chapter 
examines  the  sensing  phenomena  and  determines  a  more  precise  computational  model  than  that  in  (4.1). 
The  gains  include  less  reliance  of  calibration  procedures  in  the  reconstruction,  as  well  as  higher  quality  of 
image  reconstruction.  The  higher-order  model  is  tested  through  extensive  simulations  and  experimentally 
in  a  CAS  SI  multi-frame  testbed. 


4.2  System  Model 

The  CASSI  optical  architecture  proposed  in  [?]  is  depicted  in  Fig.  4.1.  Denote  the  spatio-spectral  power 
source  density  as  fo(x,y,\)  where  x  and  y  index  the  spatial  coordinates  and  A  indexes  the  wavelength. 
The  spatio-spectral  image  source  density  fo(x,y,X)  is  firstly  coded  by  a  coded  aperture  T(x,y).  The 
resulting  coded  field  fi(x,y,  A)  is  subsequently  dispersed  by  a  dispersive  element  before  it  impinges  onto 
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the  FPA  detector.  The  compressive  measurements  across  the  FPA  are  realized  by  the  integration  of  the 
field  f2(x,y,X)  over  the  detector’s  spectral  range  sensitivity.  The  spectral  density  at  the  output  of  the 
dispersive  element,  can  be  expressed  as, 

h(x,y,X)  =  j  J  T(x' ,y')f0{x' ,y' ,\)5(x''  -  (x  -  S(X)))S(y'  -  y)dx'dy' ,  (4.2) 

where  S(xf  —  (x  —  S(\)))5(yf  —  y)  represents  the  optical  impulse  response  of  the  system,  such  that  5(A)  = 
<a(A)(A  —  Ac)  is  the  dispersion  induced  by  the  prism  along  the  x-axis,  which  is  centered  at  the  wavelength 
Ac  and  has  a  dispersion  coefficient  <a(A).  The  resulting  intensity  image  at  the  FPA  is  the  integration  of 
the  field  ^(x,?/,  A)  over  the  detector’s  spectral  range  sensitivity  (A)  that  can  be  represented  as  g(x,y)  = 
fA  f2  (X">  Vi  X)d\.  Assuming  ideal  optical  elements,  the  energy  in  front  of  the  2D  FPA  can  be  expressed  as, 

g(x,y)=  [  T(x-  S(\),y)fo(x-  S(\),y,\)d\.  (4.3) 

Using  a  first  order  discretization  model  [?,  ?,  ?],  a  CASSI  measurement  at  the  ( m,n)th  pixel  is  given 
by,  gmn  =  f  f  p(m>  ni  y)g(xi  y)dxdy  +  wmn,  where  remn  represents  additive  noise  and  p(m,n;x,y)  = 
rect  (^  —  m,  ^  —  n)  accounts  for  the  detector  pixelation,  with  A  being  the  detector  pixel  pitch.  This  ap¬ 
proximation  however  is  coarse  leading  to  inter-pixel  blurring  in  the  detection.  The  goal  of  this  Chapter  is 
thus  to  develop  a  more  precise,  higher  order  computational  model  that  mitigates  the  inter-voxel  interfer¬ 
ence,  which  in  turn  leads  to  higher  quality  spectral  imaging.  At  the  same  time,  the  higher-order  precision 
model  allows  less  reliance  on  calibration  corrections.  To  this  end,  the  pixelation  function  p(m,n;x,y)  is 
replaced  by  defining  the  integration  limits  at  the  detector  as, 

r(ji~\~  1)A  r(m- fl)A  n 

gmn=  /  /  T(x  -  S(X),y)f0(x  -  S(X),y,X)dXdxdy,  (4.4) 

JnA  JmA  J  A 

and  the  new  discretization  model  is  then  derived  as  follows.  The  source  fo(x,y,  A)  is  discretized  as  the 
signal  fijk  where  i,  j,  and  k  are  the  discrete  indices  accounting  for  x,  y,  and  A  respectively.  A  voxel  fijk 
represents  the  intensity  concentrated  in  a  specific  spatio-spectral  region  where  x  G  [iA,  {i  +  1) A] , 
y  G  [j A,  ( j  +  1) A] ,  and  A  G  [A&,  A^+i].  Specifically,  the  ( ijk)th  voxel  is  given  by 

Afc+i  (j  +  l) A  (i+l)A 

fijk  =  /  /  /  fo(x,y,  X)dxdyd\  (4.5) 

A  k  j  A  iA 

fo(x,y,X)dxdydX  =  cijk  ■  fo(xi,yj,  Xk), 

^ ij  k 

where  represents  the  quadrature  weight,  and  x^  yj ,  and  A&  are  average  values  in  Notice  in 

(4.5)  that  the  spectral  axis  A  has  been  discretized  in  L  spectral  bands.  The  range  of  the  kth  spectral 
band  is  [A&  A^+i]  where  A^  is  the  solution  to  the  equation  given  by 

S(Afc)  -  5(A0)  =  fcA,  k  =  0, . . . ,  L  -  1,  (4.6) 

where  again  5(A)  =  <a(A)(A  —  Ac)  is  the  dispersion  of  the  prism.  Equation  (4.6)  establishes  that  the 
spectral  axis  resolution  is  determined  by  the  prism  response  5(A)  and  by  the  pitch  of  the  detector  A. 

Aperture  Imaging  /'Dispersive  Focal  Plane\ 

code  lens  I  element  ^Array  ■ 


T(x,y)  8, 

Figure  4.1:  Optical  elements  present  in  CASSI. 
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Equation  (4.5)  also  establishes  the  spatial  resolution  which  is  determined  by  the  pitch  of  the  detector  A 
that  is  assumed  to  be  equal  to  the  coded  aperture  pitch.  Using  the  ranges  of  the  L  spectral  bands  defined 
in  (4.6),  Eq.  (4.4)  can  be  expressed  as 


(n+l)A  (m+l)A 


9mn~ 


nA  rnA 


J  j  E  J  T(x-S(\),y)f0(x-S(\),y,\)d\ 

i  A  ™  A  k= 0  \  _ 


dxdy. 


(4.7) 


To  better  illustrate  the  discretization  of  the  source,  Fig.  4.2  depicts  the  physical  phenomena  described 
in  Eq.  (4.7).  Notice  that  each  cubic  voxel  when  sheared  by  the  prism  effects,  turns  into  an  oblique  voxel. 
The  oblique  voxel  is  stretched  along  the  x  axis,  such  that  when  it  is  projected  onto  the  detector  grid,  it 
impinges  onto  several  detector  pixel  elements  at  once.  Hence,  several  voxels  at  the  source  will  impact 
each  of  the  FPA  pixels. 


Dispersive  Element  Detector 


Figure  4.2:  CASSI  integration  model.  A  voxel  of  the  data  cube  is  coded  by  the  aperture  code,  sheared 
by  the  dispersive  element  with  dispersion  S( A&)  and  projected  onto  several  pixels  of  the  detector. 


The  inter-voxel  interference  is  next  examined  by  two  discretization  models  as  illustrated  in  Fig.  4.3, 
when  viewed  from  the  top  of  the  datacube.  A  linear  dispersive  function  with  slope  equal  to  one  (ds  =  1) 
is  assumed  in  the  figure.  Figure  4.3(a)  depicts  the  traditional  discretization  approach  proposed  in  [?]. 
Figure  4.3(b)  shows  the  higher  precision  discretization  model  where  the  FPA  pixel  detector  captures  en¬ 
ergy  from  several  voxels  simultaneously.  The  dotted  lines  indicate  the  spatio-spectral  region  of  the  data 
cube  integrated  in  the  (m,  n)th  detector  pixel.  It  can  be  observed  that  the  higher  precision  discretization 
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(a)  First  order  precision  model  (b)  Higher  order  precision  model 

Figure  4.3:  (Color  online)  (a)  First  order  discretization  model.  A  voxel  impinges  onto  a  single  FPA  pixel 
detector;  (b)  higher  order  discretization  model.  A  voxel  impinges  onto  three  FPA  pixels.  Notice  that  the 
light  dispersion  path  is  on  the  (A,x)  axis  (top  view). 
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of  the  dispersive  curve  leads  to  more  voxels  superimposing  in  the  formation  of  gmn. 


Figure  4.4  shows  a  zoomed  version  of  one  voxel  of  the  source  after  it  is  sheared  by  the  prism.  No¬ 
tice  that  its  energy  will  impinge  on  up  to  three  different  FPA  pixels  when  ds  =  1.  Each  voxel  at  the 
source  can  then  be  partitioned  into  three  different  regions  denoted  as  Ro,  Ri,  and  R2.  Depending  on  the 
nature  of  5(A),  a  voxel  may  affect  more  than  3  detector  elements.  Therefore,  for  a  general  dispersion 
curve,  each  of  the  integrals  in  (4.7)  can  be  rewritten  as 


(n+l)A,(m+l)A,Afc+i 


jjjT(x-S(X),y)Mx-SW,y,X)dX^y=Y.  JJJ  T(x-S( A),  y)fo(x  —  S(X),  y,  \)dxdyd\,  (4.8) 

^  A  ™  A  \_  U=0  x  _  r  o/w  ^  D 


nA,mA,Afc 


Xk{x-S(X),y}eRu 


where  d  =  ds  + 1  for  linear  dispersion,  and  d  =  max[(ra  +  l)A  —  5(A&)]  when  a  prism  exhibits  a  non-linear 
response.  Further,  let  the  discrete  version  of  the  aperture  code  T(x,  y )  be  and  using  the  representation 
in  Eq.  (4.5)  for  /o(x,  y,  A),  then 


^k  + 1 

/  // 


T(x  5(A),t/)/o(x  S(\),y,  \)dxdyd\  ^mn&;u^(ra— /c— u)n/(m— /c— u)nfc 

X k  {x  —  S(X),y}eRu 


(4.9) 


where  the  proportion  of  the  voxel  f(m-k-u)nk  contained  in  Ru  is  taken  into  account  by  the  constant  wmnku- 
Notice  that  the  subindex  k  in  w^nku  refers  to  the  spectral  interval  [A&  Afc+i].  While  the  weights  wmnku 
can  be  estimated  using  a  calibration  process,  they  can  also  be  numerically  approximated  assuming  that 
the  the  spectral  information  is  uniformly  distributed  in  the  region  delimited  by  More  specifically, 

they  are  calculated  as 


where  Ru  is  taken  in  the  respective  interval  [A&  Afc+i].  In  practice,  the  sections  Ru  can  be  calculated 
by  estimating  the  prism  response  5(A)  and  the  misalignment  between  the  coded  aperture  and  the  FPA 
detector.  Using  Eq.  (4.8)  and  (4.9),  Eq.  (4.7)  can  be  expressed  as 


w  =  E  E  ™mn*ut{m-k-u)nf(m-k-u)nk. 

fc=0  u= 0 


(4.11) 


Equation  (4.11)  can  be  written  in  matrix  form  as  g  =  Hf ,  where  the  N(N  +  L  +  d  —  1)  long  vector  g  and 
the  N2L  long  vector  f  represent  the  compressive  measurements  and  the  spectral  data  cube  respectively, 
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Figure  4.4:  A  voxel  dispersed  into  the  regions  Ro?  Ri,  and  R 2  in  each  interval  [A&  A^+i].  These  regions 
determine  the  voxel  fractions  involved  in  the  formation  of  the  9m, n  and  gm-\-i,n  detector  pixels. 
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ordered  lexicographically.  When  several  FPA  measurements  are  captured  each  one  using  a  different 
aperture  code,  the  ith  FPA  measurement  can  be  written  as 

g*  =  H,f.  (4.12) 

In  Eq.  (4.12)  each  N(N  +  L  +  d  —  1)  x  N2L  matrix  is  composed  by, 

H,  =  PT;,  (4.13) 

where  P  accounts  for  the  weights  remn/eu  and  the  dispersion  of  the  prism,  and  T \  represents  the  ith  coded 
aperture.  Here  is  a  block- diagonal  matrix  of  the  form, 

"diag(ti)  0tv2  •••  0N2 

Ow2  diag(ti)  •  •  •  0^2 

Ti  =  .  .  (4.14) 

_  0^2  0at2  •••  diag(t*)_ 

where  t i  is  the  ith  aperture  code  in  lexicographical  notation  and  0^2  is  a  N 2  x  N 2  zero-matrix.  Notice 
that  Eq.  (4.10)  can  be  written  in  matrix  form  as,  (WJf)mn  =  cjmn/cU,  for  m,  n  =  0, 1, . . . ,  N  —  1,  and  &,  u 
as  in  Eq.  (4.11).  Then,  the  matrix  P  is  given  by  P  =  Ylt=o  such  that 

®NuxN2L 

diag(Wo)  OatxTV2  **•  ®NxN2 
Otvxat2  diag(Wi)-  •  •  0NxN2  | 

Pu=  .  ....  (4.15) 

Oatxtv2  0NxN2  •  • -diag(W^_i) 

® N (d—u)  x  N2  L 

The  ensemble  of  measurements  {g^}^,  can  be  succinctly  written  as  g  =  Hf  =  PTf,  where  g  = 
[gi\  ^  g k\  ,  V  is  a  K- times  block- diagonal  matrix  of  P,  and  T  =  [Tf , . . . ,  T^]  ,  with  K  representing 
the  number  of  FPA  shots.  Figure  4.5  depicts  the  structure  of  the  matrix  H  for  the  first  order  and  the 
higher  order  discretization  models,  when  three  FPA  shots  are  used  to  capture  a  6  x  6  x  5  datacube. 

4.3  Simulations 

In  order  to  compare  the  higher-order  precision  model  with  the  traditionally  used  model,  a  hyper-spectral 
data  cube  was  experimentally  acquired  using  a  wide-band  Xenon  lamp  as  the  light  source,  and  a  visible 
monochromator.  Monochromatic  images  were  captured  every  lnm  in  the  spectral  range  {450  —  620}; 
thus  170  spectral  planes  were  acquired.  The  image  intensity  was  captured  by  a  CCD  camera  AVT 
Marlin  F033B,  with  656  x  492  pixels,  exhibiting  a  pixel  pitch  of  9.9/tm  and  using  8  bits  for  pixel  depth. 
In  addition,  a  double  Amici  prism  was  used  as  the  dispersive  element.  Its  non-linear  dispersion  curve 
shown  in  Fig.  4.10(b)  was  determined  experimentally  by  using  the  monochromator  as  the  input  of  the 
setup.  Other  elements  such  as  the  lens,  the  spectral  response  of  the  camera,  and  the  coded  apertures  are 
considered  ideal.  Deviations  from  the  ideal  characteristics  of  these  elements  are  mitigated  partially  by 
a  calibration  step.  In  this  architecture,  a  voxel  spanning  any  of  the  following  wavelength  intervals  will 
create  a  displacement  of  a  pixel  on  the  detector:  {450  —  463},  {464  —  477},  {478  —  493},  {494  —  510}, 
{511  —  530},  {531  —  556},  {557  —  586}  and  {587  —  620}  where  all  the  intervals  are  given  in  nanometres. 
Thus,  the  170  spectral  planes  of  the  datacube  will  be  clustered  into  8  bands.  Notice  that  the  intervals 
width  are  not  constant,  as  a  non-linear  dispersive  element  was  used.  Figure  4.6  shows  the  eight  spectral 
bands  of  the  datacube,  which  the  CASSI  system  aims  at  recovering. 

Experiments  use  the  multi  frame  CASSI  setup  described  in  [?],  and  simulation  algorithms  utilizes 
64  bits  as  arithmetic  precision.  Aperture  codes  entries  are  random  realizations  of  a  Bernoulli  random 
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Figure  4.5:  Structure  of  the  matrix  H  for  a  N  =  M  =  6,  L  =  5  datacube,  when  K  =  3  for  (a)  CASSI 
traditional  model  (H  G  IR180x180);  (b)  Higher  order  CASSI  model  (H  G  R180x180).  Extra  diagonal  terms 
account  for  the  inter- voxel  interference.  Notice  that  entries  in  (a)  are  either  0  or  1,  while  in  (b)  they  vary 
in  the  interval  [0,  1]. 


Figure  4.6:  Spectral  bands  used  in  the  simulations  and  their  central  wavelength. 

variable  with  parameter  p  =  0.5.  Note  that  the  proposed  model  uses  the  same  pitch  for  coded  aperture 
features  and  FPA  pixels.  In  summary,  the  hyper- spectral  test  data  cube  F  has  256  x  256  pixels  of  spatial 
resolution  and  L  =  8  spectral  bands  in  the  range  450  nm  to  620  nm.  In  order  to  simulate  the  analogous 
sensing  process,  the  compressive  measurements  are  obtained  using  the  170  spectral  planes  of  the  dat¬ 
acube,  as  shown  in  Fig.  4.3(b).  Notice  also  that  the  reconstruction  process  aims  to  recover  the  average  of 
the  spectral  information  in  the  8  mentioned  spectral  intervals.  The  calibration  weights  for  the  proposed 
model  are  approximated  using  Eq.  4.10  and  the  prism’s  response  curve.  Given  the  set  of  compressive 
measurements,  the  voxels’  weight  distribution  and  the  set  of  coded  apertures,  the  hyper- spectral  dat¬ 
acube  is  recovered  using  the  GPSR  algorithm  [?].  GPSR  exploits  the  sparse  nature  of  the  hyperspectral 
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datacube.  In  particular,  the  hyperspectral  signal  F  E  RiVxMxLj  or  its  vector  representation  f  E  RAr  M  L, 
are  assumed  to  be  K— sparse  on  some  basis  such  that  f  =  ^3 £>#,  where  6  are  the  coefficients  of  the 
sparse  representation.  Hence,  f  can  be  approximated  by  a  linear  combination  of  K  vectors  from  with 
K  <C  ( N.M.L ).  Specifically,  this  algorithm  estimates  a  hyperspectral  datacube  f  by  solving  the  optimiza¬ 
tion  problem,  f  =  ^^arguiin^, \\g  —  VT'b 3D^7 1| i  +  Tll^/||i}>  where  r  >  0  is  a  regularization  parameter 
that  balances  the  conflicting  tasks  of  minimizing  the  least  square  of  the  residuals,  while  at  the  same  time, 
it  seeks  for  a  sparse  solution  The  basis  representation  is  set  as  the  Kronecker  product  of  three  basis 
^3  d  =  ^ri(8)^r2^^r3,  where  the  combination  ®  is  the  2D- Wavelet  Symmlet  8  basis  and  ^3  is  the 
Discrete  Cosine  basis.  The  reconstructions  are  performed  using  the  new  model  in  Eq.  (4.11),  and  the 
traditional  model  in  Eq.  (4.1)  with  its  respective  calibration  process  described  in  [?].  The  regularization 
parameters  needed  in  the  compressive  sensing  reconstruction  algorithm  are  carefully  selected  such  that 
each  simulation  uses  the  best  selectable  parameter.  Figure  4.7  depicts  the  reconstructed  spectral  bands 
(zoomed  area)  when  6  shots  are  captured  for  the  model  in  (4.1).  Figure  4.8  illustrates  the  reconstruction 
of  the  same  spectral  bands  (zoomed  area)  when  the  same  number  of  shots  are  used  in  the  new  model.  It 
can  be  observed  that  the  new  model  recovers  the  spectral  information  with  higher  accuracy. 
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Figure  4.7:  Reconstruction  using  the  traditional  CASSI  model  and  the  corresponding  attained  PSNR. 
The  average  PSNR  across  the  8  bands  is  22.3  dB. 


Figure  4.9  shows  the  PSNR  of  the  reconstructions  for  the  two  models  as  function  of  the  measurement 
shots.  The  gain  achieved  by  the  new  model  is  quantitatively  noticeable  by  averaging  the  PSNR  of  the 
recovered  datacubes.  This  improvement  approaches  to  4  dB  when  more  than  two  FPA  shots  are  used. 

4.4  Experimental  results 

The  testbed  shown  in  Fig.  4.10(a)  is  used  to  implement  the  CASSI  system  and  to  verify  the  simulation 
results  [?].  It  is  formed  by  two  subsystems:  the  first  composed  by  the  illuminated  target,  the  objective 
lens  and  the  DMD;  the  second  by  the  imaging  lenses,  the  band  pass  filter,  the  dispersive  element,  and 
the  CCD  camera.  The  target  is  illuminated  by  a  white  light  source  and  its  reflected  light  is  captured 
by  the  objective  lens  which  focuses  the  light  onto  the  DMD  plane,  which  plays  the  role  of  the  coded 
aperture.  Afterwards,  the  reflected  light  from  the  DMD  is  focused  by  the  imaging  lenses  into  the  prism 
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Figure  4.8:  Reconstruction  using  the  higher  order  CASSI  model  and  the  corresponding  attained  PSNR. 
The  average  PSNR  across  the  8  bands  is  26.85  dB 


Figure  4.9:  Averaged  PSNR  of  the  reconstructed  datacubes  as  function  of  the  number  of  FPA  shots.  The 
traditional  and  the  higher  order  precision  models  are  compared. 


imaging  plane  that  disperses  the  filtered  light  onto  the  CCD  camera  which  integrates  the  underlying  3D 
hyperspectral  image  in  the  2D  FPA. 

The  testbed  setup  is  characterized  in  order  to  reduce  the  impact  of  non-linearities,  non- uniformities,  and 
external  noise  artifacts.  This  process  is  realized  as  follows:  (a)  The  light  source  intensity  distribution  and 
the  FPA  spectral  sensitivity  are  characterized  by  experimentally  analyzing  their  spectral  responses  using 
a  USB2000+VIS-NIR  Ocean  Optics  spectrometer  with  a  known  spectral  response.  These  non-uniform 
spectral  response  curves  are  taken  into  account  to  reduce  their  impact  in  the  measurement  shots;  (b)  for 
each  one  of  the  170  captured  spectral  planes,  10  FPA  measurements  are  captured  and  averaged  to  reduce 
the  impact  of  shot  and  readout  noise;  (c)  the  CCD  exposure  time  is  setted  to  100  microseconds,  in  order 
to  improve  the  signal-to-noise-ratio  of  the  aperture  code  at  each  wavelength;  (d)  the  dispersive  element  is 
characterized  in  order  to  take  into  account  its  non-linear  response  curve  and  the  resultant  bandwidth  of 
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(a)  CASSI  testbed  setup 


(b)  Non-linear  response  of  the  Amici  prism 


Figure  4.10:  (a)  The  CASSI  testbed  setup  and  its  six  optical  elements:  objective  lens,  DMD,  imaging 
lenses,  band-pass  filter,  prism  and  CCD;  (b)  non-linear  dispersion  response  of  the  Amici  prism  between 
{450  -  620}nm. 


each  spectral  band.  After  characterization  of  the  testbed  and  in  order  to  observe  the  oblique  voxel  effect 
impinging  into  the  FPA  as  explained  in  Fig.  4.3,  a  measurement  shot  is  captured  using  monochromatic 
light  at  502  nm  as  the  input  of  the  system.  The  resultant  measurement  is  depicted  in  Fig.  4.11. 


Coded  aperture  FPA  measurement 


Figure  4.11:  FPA  measurement  at  502  nm.  The  coded  aperture  (upper- left)  is  used  in  order  to  isolate 
the  effect  of  a  single  voxel  impinging  onto  the  FPA  (upper  right).  A  zoomed  version  of  a  single  FPA 
pixel  shows  the  measured  intensity  taken  into  account  for  each  of  the  discretization  models.  The  energy 
classified  as  noise  and  blur  by  the  first  order  and  the  higher  order  models,  is  shown. 

This  measurement  is  taken  using  a  test  coded  aperture  with  enough  space  (3  ‘off’  features)  between 
each  ‘on’  feature,  thus,  allowing  the  isolation  of  the  effect  of  a  single  voxel  impinging  onto  the  FPA.  Then, 
a  zoomed  version  of  a  single  FPA  pixel  is  analyzed.  Firstly,  it  can  be  confirmed  that  energy  belonging 
to  a  single  datacube  voxel  impinges  principally  in  a  single  FPA  pixel  (m,  n)t/l,  and  a  smaller  portion  is 
projected  into  its  neighbors  (m  —  1  ,n)th  and  (m  +  1  ,n)th.  Second,  the  first  order  model  accounts  only 
for  the  energy  impinging  on  the  principal  pixel,  discarding  the  energy  around  it,  or  classifying  it  as  noise 
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Figure  4.12:  Objects  in  scene  used  in  the  experimental  comparison 
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Table  4.1:  Weights  Ri  and  their  distribution  across  spectral  bands 


or  blur.  The  energy  discarded  by  the  first  order  model,  is  taken  into  account  in  the  higher  order  model 
by  the  weights  remn/cu.  Notice  that  the  energy  considered  as  noise  and  blur  by  the  traditional  CASSI  is 
leveraged  by  the  use  of  a  calibration  cube  at  the  reconstruction  stage,  while  the  proposed  model  finds 
the  weights  distribution  in  an  off-line  process.  Then,  the  proposed  higher-precision  computational  model 
becomes  more  suitable  for  reconfigurable  multi-shot  CASSI  where  multiple  coded  apertures  are  used  se¬ 
quentially. 

For  experimental  purposes,  the  objective  scene  used  is  depicted  in  Fig.  4.12.  The  coded  apertures 
are  realizations  of  a  Bernoulli  random  variable  with  p  =  0.5,  realized  by  a  DMD  as  explained  in  [?] 
exhibiting  128  x  128  pixels.  The  dispersive  element  is  an  Amici  prism  exhibiting  the  non-linear  response 
shown  in  Fig.  4.10(b).  To  match  with  the  pitch  of  the  coded  aperture  features  and  accounting  for  the 
dispersion  process,  128  x  136  pixels  of  the  CCD  are  required.  The  weight  distribution  extraction  is  per¬ 
formed  using  Eq.  (4.10).  As  a  result,  three  128  x  128  x  8  weight  datacubes  were  obtained,  each  one 
accounting  for  the  regions  i?o,  Ri  and  R 2  as  described  in  Fig.  4.4.  By  averaging  each  weight  datacube 
per  band,  a  succinctly  version  can  be  shown  in  Table  4.1.  Notice  that  when  misalignments  between  the 
coded  aperture  and  the  FPA  occur,  the  weights  distributions  may  vary  along  the  regions. 

The  GPSR  algorithm  is  employed  in  the  reconstruction  of  the  underlying  hyperspectral  scene,  with  pa¬ 
rameters  as  described  in  the  simulations  section  [?].  Figure  4.13  depicts  the  8  reconstructed  spectral 
bands  for  6  shots  when  the  models  in  Eq.  (4.1)  and  Eq.  (4.11)  are  used.  Here,  the  higher  quality  recon¬ 
struction  obtained  when  the  proposed  model  is  used  in  the  simulations  section  is  confirmed.  Notice  that 
the  test  object  intensity  spans  principally  along  the  last  four  bands,  and  the  reconstruction  quality  of  the 
proposed  model  overcomes  the  one  from  the  traditional  CASSI  model.  In  particular,  the  improvement 
can  be  clearly  noticed  in  the  fifth  band  (524  nm),  where  the  higher-order  CASSI  estimates  a  better  shape 
of  the  face,  compared  with  the  same  band  from  the  traditional  CASSI.  Furthermore,  the  improved  results 
can  also  be  noticed  in  the  spectral  signatures  of  two  particular  points  (PI  and  P2  in  Fig.  4.12)  depicted 
in  Fig.  4.14.  The  resulting  reconstructed  datacubes  curves  are  compared  against  their  respective  ground 
truth  curves  measured  by  the  use  of  the  commercial  spectrometer. 

It  is  important  to  point  out  here,  that  the  simulations  setup  differs  from  the  experimental  setup  in  the 
following  aspects:  The  former  performs  the  CASSI  and  higher-order  CASSI  FPA  measurements  starting 
with  a  hyperspectral  datacube  captured  off-line  and  taken  as  the  ground  truth.  The  coded  apertures  as 
well  as  the  non-linear  prism  dispersion  curve  are  simulated,  and  the  weights  distribution  (wmnku)  given 
by  Eq.  (4.10)  are  then  synthetically  obtained.  Noise,  as  well  as  blur  and  misalignment  between  the 
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coded  aperture  and  the  CCD  were  not  added  to  the  model.  It  can  be  assumed  that  this  is  the  ideal  case 
scenario.  The  latter  uses  the  experimental  testbed  depicted  in  Fig.  4.10(a)  to  capture  the  FPA  mea¬ 
surement  shots.  The  same  fluorescent  white  illumination  source  was  employed  for  both  cases.  However, 
in  the  experimental  results,  the  FPA  measurements  are  contaminated  by  optical  aberrations,  as  well  as 
noise  and  misalignment  between  the  CCD  and  the  coded  aperture.  Consequently,  the  weights  distribu¬ 
tion  (wjnnku)  presented  in  Table.  4.1,  is  experimentally  obtained  from  the  non-linear  prism  dispersion 
curve  depicted  in  Fig.  4.10(b),  which  is  characterized  by  the  use  of  a  monochromatic  light  source  as  the 
input  of  the  testbed  ranging  between  450nm  and  650nm.  The  reconstruction  algorithm  utilized  (GPSR) 
is  equal  for  both  simulations  and  experimental  results,  but  differs  for  each  compared  model  (CAS SI  vs 
Higher-Order  CASSI).  The  objects  used  as  targets  in  the  simulation  section  differs  from  the  one  used 
in  the  experimental  setup.  Due  to  field  of  view  restrictions  in  the  optical  instruments,  a  smaller  but 
spatially  richer  scene  was  selected  for  the  latter. 

4.5  Conclusions 

A  higher  order  precision  discretization  model  for  coded  aperture-based  spectral  imaging  systems  has  been 
developed.  This  model  accounts  for  the  inter-voxel  projections  onto  each  pixel  detector  which  is  disre¬ 
garded  by  the  first  order  discretization  model.  This,  in  turn,  allows  for  the  reconstruction  of  hyperspectral 
signals  with  higher  PSNR.  Simulations  achieve  a  4  dB  improvement,  while  testbed  experiments  visually 
confirm  the  simulations  results.  The  proposed  model  is  less-dependent  on  time-demanding  calibration 
processes,  thus  leading  to  multiple- frame  CASSI  systems  to  be  more  suitable  for  real  applications. 
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(a)  CASSI 


(b)  Higher-order  CASSI 

Figure  4.13:  Reconstruction  of  the  8  spectral  bands  using  (a)  the  traditional  CASSI  model,  and  (b)  the 
proposed  higher-order  CASSI  model. 
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(a)  Red  point  (PI  in  Fig.  12)  (b)  Yellow  point  (P2  in  Fig.  12) 

Figure  4.14:  Spectral  signatures  comparison  from  given  points  in  Fig.  12. 
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Chapter  5 


Compressive  Hyperspectral  Imaging 
Testbed 


In  this  chapter,  we  present  a  dispersive  optical  element-based  approach  for  building  a  CS-HSI  system. 
This  system  has  a  compact  optical/mechanical  architecture,  along  with  a  higher  efficiency  in  utilizing 
the  DMD  reflected  imaging  irradiance  and  therefore,  spectral  images  can  be  concurrently  reconstructed 
from  this  system.  Additionally,  the  acquisition  of  spatial/spectral  image  cubes  using  the  CS-HSI  system 
is  free  from  mechanical  or  temporal  scanning  processes  and  thus  this  acquisition  process  virtually  does 
not  cause  any  time  penalty. 

5.1  Coded  Aperture  Snapshot  Spectral  Imaging  System 

Conventional  spectral  imaging  systems  usually  rely  on  some  mechanical  or  temporal  scanning  processes 
to  acquire  the  complete  spatial/spectral  information  of  an  imaging  scene  [49].  The  time  penalty  caused 
by  the  scanning  process  unavoidably  undermines  the  performance  of  spectral  imaging  systems  in  low  light 
or  high  speed  imaging  applications.  To  circumvent  the  scanning  process,  the  so-called  Coded  Aperture  S- 
napshot  Spectral  Imaging  (CASSI)  system  was  developed  [50-53].  Figure  5.1  shows  the  schematic  drawing 
of  the  CASSI  system.  CASSI  systems  usually  use  photomasks  to  implement  CS  measurement  patterns. 


Figure  5.1:  Schematic  drawing  of  the  CASSI  system. 

Binary  random  patterns  are  widely  used  as  CS  measurement  patterns  in  such  systems.  The  photomask 
is  installed  on  the  image  plane  of  an  imaging  lens  such  that  intensities  of  optical  images  formed  on  the 
photomask  can  be  modulated  by  the  CS  measurement  pattern  coated  on  it.  In  this  process,  the  continu¬ 
ous  spatial  information  of  the  optical  image  is  pixelated  into  an  array  of  square  pixels  by  the  photomask. 
In  the  following  discussion,  we  use  f(m ',  n',  k)  to  represent  the  spatial/spectral  information  of  the  optical 
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image  formed  on  the  photomask,  and  use  t(m',n')  to  represent  the  CS  measurement  pattern.  The  coor¬ 
dinates  m'  and  n'  are  used  to  locate  a  pixel  on  the  photomask  ( m '  =  1,2,...,  M'  and  n'  =  1,  2, . . . ,  TV7). 
The  coordinate  k  is  used  to  describe  the  location  of  a  spectral  image  in  the  wavelength  (A)  direction  of  the 
image  cube  {k  =  1,2,.. . ,  if).  The  spatial/spectral  information  of  the  original  optical  image  after  being 
modulated  by  the  CS  measurement  pattern  can  be  represented  as:  o(m',n',fc)  =  f(m',n',k)t(m',n').  A 
relay  lens  is  used  to  transfer  the  modulated  image  intensity  from  the  photomask  onto  a  CCD  camera. 
In  between  the  exit  aperture  of  the  relay  lens  and  the  CCD  camera,  an  optical  dispersive  element  (such 
as  a  prism)  is  used  to  spatially  separate  the  spectral  components  of  the  modulated  image  intensity.  In 
other  words,  a  spectrally  dispersed  image  of  the  modulated  image  intensity  is  formed  on  the  CCD  cam¬ 
era,  which  represents  a  set  of  CS  measurement  results  for  the  original  image  cube.  We  use  g(m,n )  to 
represent  that  spectrally  dispersed  image,  where  (m,  n)  is  the  coordinate  system  used  to  locate  a  pixel  on 
the  optical  plane  where  the  CCD  camera  is  installed.  If  the  magnification  of  the  relay  lens/double- Amici 
prism  structure  is  A,  we  have  m  =  Am'  and  n  =  An' .  In  following  discussions,  we  assume  A  =  1. 

In  [53],  optical  operations  realized  by  the  CASSI  system  were  described  as  punch,  shear,  and  smash. 
Figure  5.2  provides  an  illustrative  explanation  for  those  optical  operations.  In  Fig.  5.2,  we  only  show 
three  spectral  channels  in  the  original  image  cube,  including  a  red  channel,  a  green  channel,  and  a  blue 
channel.  The  spatial  shift  of  different  spectral  channels  caused  by  the  prism  happens  only  in  the  x  direc¬ 
tion  on  the  CCD  camera  plane  and  we  call  that  direction  the  dispersion  direction  of  the  CASSI  system. 
Also  in  Fig.  5.2,  the  green  spectral  channel  is  spatially  shifted  from  the  red  channel  by  1  CCD  pixel 
to  the  right  in  the  dispersion  direction,  and  the  blue  channel  is  shifted  from  the  red  channel  by  2  CCD 
pixels  to  the  right. 
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CS  measurement 
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Detector  pixel  lattice 
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3-D  spatial/ spectral  Punch 

image  cube 


Figure  5.2:  Optical  operations  realized  by  the  CASSI  system. 

The  first  optical  operation  realized  by  the  CASSI  system  is  called  punch,  wherein  closed  apertures  in 
the  photomask  punch  the  spatial/spectral  information  away  from  the  original  image  cube  in  corresponding 
pixel  locations.  The  second  operation  is  called  shear,  wherein  different  spectral  components  of  the  punched 
image  cube  are  spatially  separated  from  each  other  along  the  dispersion  direction,  due  to  the  dispersion 
performance  of  the  prism.  The  last  optical  operation  is  called  smash,  wherein  the  continuous  information 
of  the  sheared  image  cube  is  digitally  sampled  and  collected  by  a  CCD  camera.  Those  optical  operations 
were  modeled  as  a  system  transfer  function  H  in  [51-53].  If  we  use  g  to  represent  the  spectrally  dispersed 
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image  captured  by  the  camera,  the  image  formation  process  of  the  CAS  SI  system  can  be  modeled  as: 

9  =  Hf.  (5.1) 


5.2  Transfer  Matrix  of  the  CASSI  System 

In  Fig.  5.2,  we  note  that  the  intensity  of  the  CS  measurement  result  g  at  the  pixel  location  (m,  n)  equals 
to  the  intensity  of  the  first  spectral  channel  (k  =  1)  of  the  modulated  image  intensity  at  the  pixel  location 
(m,  n)  plus  the  intensity  of  the  second  spectral  channel  (k  =  2)  at  the  pixel  location  (m,  n  —  1),  plus  the 
intensity  of  the  third  spectral  channel  (k  =  3)  at  the  pixel  location  (m,  n  —  2),  and  so  on.  Following  this 
logic,  for  an  M  x  TV  x  K  image  cube,  the  CS  measurement  result  generated  from  the  CASSI  system  can 
be  represented  as: 


Ylk= l  °(m > n  ~  k  l,k),  if  1  <  n  <  K  —  1 

g(m ,  n)  =  \  Yk=%  °(m > n  ~  k  +  1,  k),  if  K  <  n  <  TV 


(5.2) 


=n  °(mi  n~  k  +  l,k),  ifN  +  l<n<N  +  K  —  1 


where  o(m,  n  —  k  +  1,  k)  =  /(m,  n  —  k  +  1  )t(m,  n  —  k  +  1).  To  determine  the  system  transfer  matrix  H ,  we 
transform  the  2D  and  3D  data  arrays  in  Eq.  5.2,  which  are  g  and  o,  respectively,  into  ID  vector  forms  by 
concatenating  their  column  vectors.  In  their  ID  vector  forms,  g  has  M  x  (TV  +  K  —  1)  data  entries  and 
/  has  M  x  N  x  K  data  entries.  Thus,  the  transfer  matrix  H  has  a  pixel  dimension  of  M  x  (TV  +  K  —  1) 
by  M  x  A  x  if,  which  is  a  rectangle  with  an  aspect  ratio  smaller  than  1.  We  define  the  aspect  ratio  of 
the  transfer  matrix  to  be  the  down-sampling  ratio  of  the  CASSI  system,  which  is  denoted  as  Rcassi  in 
this  work,  and  we  have: 

M  x  (TV  +  K  —  1)  N  +  K  —  1 


Rcassi  = 


(5.3) 


M xN xK  NxK 

Apparently,  the  down-sampling  ratio  decreases  with  the  increment  of  the  number  of  spectral  channels 
in  the  image  cube.  In  their  ID  vector  representations,  the  pixel  location  (m,  n)  in  g  becomes  the  ((n  — 
1  )M  +  m)th  data  entry  and  the  pixel  location  (m,  n,  k)  in  /  becomes  the  ((k  —  1)MN  +  (n  —  1  )M  +  m)th 
data  entry.  The  pixel  location  (m,  n  —  k  +  1,  k)  in  /  becomes  the  ((fc  —  1)MTV  +  (n  —  k)M  +  m)th  data 
entry.  If  we  use  the  ID  coordinate  ((n  —  1  )M  +  m)  to  replace  the  2D  coordinate  (m,  n)  in  Eq.  5.2,  use 
((k  —  1)MN  +  (n  —  1  )M  +  m)  to  replace  (m,  n,  fc)  and  let  z  m  (n  —  1)M  +  m,  which  represents  the  pixel 
index  in  the  ID  form  of  the  CS  measurement  result  g,  the  center  term  of  Eq.  5.2  can  be  modified  as: 

K 


g(z)  =  o((k  —  1)MN  +  (n  —  k)M  +  m) 


k=  1 
K 


=  o((fc  —  1)MTV  +  nM  —  kM  +  m) 

k=l 

K 

=  ^  o((fe  -  1)MTV  +  nM  -M  +  M  -kM  +  m) 

k=  1 
K 

=  °((k  -  l)MN  +  z-(k-  l)Af ), 


(5.4) 


k= 1 


which  is  valid  when  (if  —  1)M  +  1  <  2  <  NM.  Similarly,  the  top  term  in  Eq.  5.2  can  be  modified  as: 

r*i 

g(z)  =  Y  °((k  ~  l)NM  +  ^  -  (*  -  1)M),  (5.5) 


k= 1 


which  is  valid  when  1  <  z  <  ( K  —  1  )M.  \jf  \  represents  the  ceiling  function  of  The  bottom  term  in 
Eq.  5.2  can  be  modified  as: 


K 


g(z)  =  °((fc  -  l)NM  +  z  —  (k  —  1  )M), 

r^i 


(5.6) 
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which  is  valid  when  MN  +  1  <  z  <  M{N  +  K  —  1). 

The  explicit  expression  of  the  transfer  matrix  H  can  be  derived  from  Eq.  5.4  through  Eq.  5.6.  Here 
we  use  a  simple  example  to  demonstrate  the  derivation  process,  where  we  consider  a  small  image  cube 
which  has  M  =  N  =  3  and  K  =  3.  The  CS  measurement  pattern  t  has  a  pixel  dimension  of  3  x  3,  and 
the  CS  measurement  result  g  has  a  pixel  dimension  of  3  x  5.  In  their  ID  vector  forms,  t  and  g  have  9  and 
15  data  elements,  respectively.  When,  1  <  z  <  6,  Eq.  5.5  should  be  used.  When  z  =  1,  \jf  \  =  1  and  we 
have  g(  1)  =  J2k=i  =  °((&  —  1  )NM  +  1  —  {k  —  1  )M)  =  o(l)  =  £(1)/(1).  Similarly,  we  have  g(2)  =  £(2)/(2) 
and  g( 3)  =  £(3)/(3).  In  other  words,  data  entries  in  pixel  locations  of  (1,1),  (2,2),  and  (3,3)  in  the 
transfer  matrix  H  are  t(  1),  t( 2),  and  t( 3),  respectively.  The  last  are  the  first  three  data  entries  in  the 
ID  vector  form  of  the  CS  measurement  pattern.  When  z  =  4,  also  from  Eq.  5.5,  |~ =  2  and  we  have 
g( 4)  =  Y^k= i  °((k  ~  1  )NM  +  4  —  (fc  —  1)M)  =  o(4)  +  o(10)  =  £(4)/(4)  +  £(10)/(10).  £(10)  represents  the 
intensity  modulation  implemented  on  the  first  pixel  location  of  the  second  spectral  channel  in  the  original 
image  cube.  Since  in  the  punch  process,  the  same  CS  measurement  pattern  was  used  to  implement  the 
intensity  modulation  on  all  the  spectral  channels  of  the  original  image  cube,  the  intensity  modulation 
implemented  on  the  first  pixel  location  in  the  second  spectral  channel  equals  to  the  intensity  modulation 
implemented  on  the  first  pixel  location  in  the  first  spectral  channel  and  thus  £(10)  =  £(1).  Similarly,  we 
have  g( 5)  =  £(5)/(5)  +  £(2)/(ll),  and  g(6)  =  £(6)/(6)  +  £(3)/(12).  In  other  words,  data  entries  in  pixel 
locations  of  (4,4),  (4,10),  (5,5),  (5,11),  (6,6),  (6,12)  in  the  transfer  matrix  are  £(4), £(1), £(5), £(2), £(6), 
and  t{ 3)  respectively.  Other  data  entries  in  the  transfer  matrix  can  be  derived  in  a  similar  manner  and 
the  resulted  transfer  matrix  H  is  shown  in  Fig.  5.3.  The  values  of  those  empty  entries  in  this  transfer 
matrix  are  all  0. 
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Figure  5.3:  Transfer  matrix  of  the  CAS  SI  system  for  an  M  =  N  =  K  =  3  image  cube. 

For  an  M  x  N  x  K  image  cube,  the  transfer  matrix  looks  like  this: 

We  can  see  that  from  row  (K  —  1)M  +  1  to  row  NM ,  each  row  in  the  transfer  matrix  contains  a  random 
data  sequence  that  has  K  data  entries.  For  rows  1  to  (K  —  1  )M  and  rows  MN -\- 1  to  M(N  +  K  —  1),  each 
row  in  the  transfer  matrix  contains  a  random  data  sequence  that  has  \jf\  and  K  —  \jf\  data  entries, 
respectively,  where  z  represents  the  row  number  of  the  H  matrix,  which  is  identical  to  the  row  number  of 
the  CS  measurement  result  g.  As  such,  each  row  in  the  transfer  matrix  H  implements  a  random  intensity 
modulation  (random  sampling)  on  the  original  image  cube  /,  and  each  pixel  in  the  spectrally  dispersed 
image  g  represents  a  CS  measurement  result  for  /.  In  CASSI  systems,  random  sampling  processes  are 
implemented  on  the  3D  spatial/spectral  image  cube  of  optical  images.  We  utilized  a  Two  Step  Iterative 
Shrinkage/Thresholding  (TwIST)  algorithm  to  solve  the  image  reconstruction  problem,  which  considers 
the  following  Total  Variation  (TV)-regularized  minimization  problem  [51-54]: 

argminfeRN  \\g  -  Hf\\\  +  rTV(f),  (5.7) 

where  r  is  a  regularization  parameter  for  the  Total  Variation  regularizer  TV(f),  which  is  defined  as: 

TV  {f)  =  EE  V {Jim  +  1,  n,  k)  -  J{m,  n,  fc))2  +  (/(m,  n  +  l,k)~  f(m,  n,  A:))2.  (5.8) 

k  m,n 
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Figure  5.4:  Transfer  matrix  for  an  M  x  N  x  K  image  cube. 


5.3  Simulations 

We  simulated  the  punch,  shear,  and  smash  operations  of  the  CAS  SI  system  in  the  simulation  setting  (we 
used  Matlab  as  the  simulation  tool)  and  used  the  TwIST  algorithm  to  reconstruct  the  original  image 
cube  based  on  the  simulated  CS  measurement  results.  In  this  simulation,  a  256  x  256  x  24  data  cube  was 
used  as  the  imaging  target,  which  is  shown  in  Fig.  5.5(a). 


Figure  5.5:  (a)  Spatial/spectral  image  cube  of  the  simulation  target,  (b)  RGB  image  of  the  simulation 
target,  (c)  Optical  image  of  the  simulation  target. 


Figures  5.5(b)  and  (c)  show  an  RGB  image  and  an  optical  image  of  the  simulation  target.  The  RGB 
image  was  generated  by  using  the  454  nm,  550  nm,  and  621  nm  spectral  channels  to  represent  the  blue, 
green,  and  red  colors  in  the  RGB  data  format.  In  the  RGB  image,  we  can  see  that  this  simulation  target 
is  composed  of  some  red  and  green  tomatoes.  The  optical  image  is  the  intensity  combination  of  all  the 
24  spectral  channels  of  the  image  cube.  In  Fig.  5.6,  we  show  the  normalized  intensity  variations  of  two 
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pixel  locations  in  the  image  cube,  which  are  labeled  as  the  red  spot  (y=36,  x=206)  and  the  green  spot 
(y=222,  x=122)  in  Fig.  5.5(b).  The  red  spot  is  located  on  a  red  tomato  and  the  green  spot  is  located  on 
a  green  tomato.  Those  intensity  variations  can  be  considered  as  the  spectral  curves  of  the  imaging  target 
at  those  two  pixel  locations.  For  the  red  spot,  the  spectral  curve  has  an  intensity  peak  at  the  spectral 
channel  20,  which  represents  a  wavelength  of  607  nm.  For  the  green  spot,  there  is  an  intensity  peak  at 
the  spectral  channel  14,  which  represents  a  wavelength  of  538  nm. 


Figure  5.6:  Spectral  curves  measured  at  the  red  spot  and  the  green  spot  on  the  imaging  target. 

A  256  x  256  binary  random  pattern  was  used  to  implement  the  CS  measurement  process,  which  is 
shown  in  Fig.  5.7(a).  In  this  CS  measurement  pattern,  half  of  its  pixels  are  set  in  the  on  condition  and 
thus  we  say  this  pattern  has  an  on  pixel  percentage  of  50%.  The  CS  measurement  result  generated  with 
this  random  pattern  is  shown  in  Fig.  4.7(b),  which  was  fed  into  the  TwIST  algorithm  for  reconstructing 
the  original  image  cube.  In  this  simulation,  we  assume  the  spectral  dispersion  caused  by  the  prism 
happens  only  in  the  horizontal  direction  on  the  CCD  camera  plane  and  adjacent  spectral  channels  are 
spatially  separated  by  1  pixel.  Therefore,  the  CS  measurement  result  has  a  pixel  dimension  of  256  x  279. 
In  the  TwIST  algorithm,  we  set  the  regularization  parameter  r  to  be  0.1  and  the  number  of  optimization 
iterations  to  be  100.  The  image  reconstruction  took  about  650  seconds  using  a  DELL  desktop  (model: 
precision  690)  installed  with  Windows  XP  32  bit  operation  system  and  4GB  memory. 


Figure  5.7:  (a)  A  256  x  256  CS  measurement  pattern,  (b)  Simulated  CS  measurement  result  for  the 
image  cube  shown  in  Fig.  5.5(a). 

Figure  5.8(a)  shows  the  reconstructed  image  cube.  Figures  5.8(b)  and  (c)  show  the  RGB  image  and 
the  optical  images  of  the  reconstructed  image  cube. 
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Figure  5.8:  (a)  Reconstructed  image  cube,  (b)  RGB  image  of  the  reconstructed  image  cube,  (c)  Optical 
image  of  the  reconstructed  image  cube. 


In  Fig.  5.8(a),  we  can  see  that  the  red  tomatoes  in  the  imaging  target  were  reconstructed  with  strong 
intensities  in  red  spectral  channels  (595  nm  to  637  nm)  and  green  tomatoes  were  reconstructed  with 
strong  intensities  in  green  spectral  channels  (529  nm  to  582  nm).  The  PSNR  value  of  the  reconstructed 
optical  image  is  17.36  dB.  In  Fig.  5.9,  we  show  normalized  intensity  variations  of  the  reconstructed  image 
cube  at  the  red  spot  and  the  green  spot  (solid  lines).  Those  intensity  variations  can  be  considered  as 
the  reconstructed  spectral  curves  at  those  two  pixel  locations.  Also  in  Fig.  5.9,  spectral  curves  of  the 
original  image  cube  are  shown  in  dotted  lines.  We  can  see  that  the  reconstructed  spectral  curves  have 
pretty  good  matches  to  the  original  ones.  The  Mean  Square  Errors  (MSEs)  of  the  deviations  from  the 
reconstructed  spectral  curves  to  the  original  spectral  curves  are  0.017%  and  0.027%  for  the  red  and  green 
spots,  respectively. 

5.4  DMD-based  CASSI  System 

Conventional  CASSI  systems  use  photomasks  to  implement  CS  measurement  patterns.  In  such  system- 
s,  the  pattern  replacement  process  is  a  very  time-consuming  process.  A  new  photomask  needs  to  be 
lithographically  fabricated  and  the  entire  system  needs  to  be  re-aligned  after  the  installation  of  the  new 
photomask.  To  expedite  the  pattern  replacement  process,  we  propose  using  a  DMD  to  implement  CS 
measurement  patterns.  As  such,  new  patterns  can  be  implemented  without  optical  and/or  mechanical 
modifications  to  the  system.  More  importantly,  the  DMD  provides  a  virtually  unlimited  selection  pool 
for  CS  measurement  patterns,  such  that  patterns  of  different  designs  can  be  grouped  to  realize  various 
Multi-Shot  CS  (MS-CS)  measurement  processes,  which,  as  will  be  demonstrated  in  the  following  section- 
s,  can  effectively  enhance  the  quality  of  the  reconstructed  image  cubes  [52,55].  Figure  5.10  shows  the 
schematic  drawing  of  the  DMD-based  CASSI  (DMD-SSI)  system. 

As  can  be  seen  in  Fig.  5.10,  the  DMD-SSI  system  can  be  considered  as  a  two-  arm  system,  including  an 
imaging  arm  and  a  relay-dispersion  arm.  In  the  imaging  arm,  a  DMD  is  installed  on  the  image  plane  of  an 
imaging  lens,  which  imposes  intensity  modulations  on  optical  images  formed  on  it.  In  the  relay-dispersion 
arm,  image  intensities  modulated  by  DMD  patterns  are  collected  into  a  relay  lens/double- Amici  prism 
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Figure  5.9:  Spectral  curves  measured  at  the  red  spot  and  the  green  spot  in  the  reconstructed  image  cube 
(solid  lines)  and  spectral  curves  measured  at  those  two  locations  in  the  original  image  cube  (dotted  lines) . 


structure.  On  the  image  plane  of  that  relay  lens/double- Amici  prism  structure,  where  a  CCD  camera  is 
installed,  spectrally  dispersed  images  of  the  modulated  image  intensities  are  formed.  Those  spectrally 
dispersed  images  represent  different  sets  of  CS  measurement  results  for  the  original  image  cube.  Each  on 
pixel  on  the  DMD  mirror  plane  results  in  a  dispersion  band  on  the  CCD  camera  and  CS  measurement 
results  are  essentially  intensity  combinations  of  dispersion  bands  at  different  spatial  locations  on  the 
CCD  camera.  The  construction  of  the  imaging  arm  of  the  DMD-SSI  system  is  straightforward,  which 
only  involves  using  an  imaging  lens  to  form  optical  images  on  the  DMD  mirror  plane.  The  construction 
of  the  relay-dispersion  arm  of  the  DMD-SSI  system  is  more  complicated  than  the  imaging  arm.  Figure 
5.11(a)  shows  the  Zemax  ray  tracing  model  of  the  relay-dispersion  arm  installed  in  our  DMD-SSI  system. 
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Figure  5.11:  (a)  Zemax  ray  tracing  model  for  the  relay-dispersion  arm  of  the  DMD-SSI  system,  (b) 
Optical  design  of  the  double-Amici  prism,  (c)  Manufactured  double-Amici  prism  and  its  mechanical 
holder. 

In  this  setup,  we  used  a  pair  of  visible  achromatic  lenses  (Thorlabs,  stock  number:  AC254-100-A- 
ML)  to  realize  the  functionality  of  optical  relay.  The  double-  Amici  prism  was  installed  in  between  the 
exit  aperture  of  the  relay  lens  and  the  CCD  camera.  Figure  5.11(b)  shows  the  optical  design  of  the 
double-Amici  prism,  which  is  composed  of  three  wedge  prisms  including  a  high  dispersion  wedge  prism 
sandwiched  in  between  two  low  dispersion  ones.  The  angles  and  thicknesses  of  the  three  wedge  prisms 
were  optimized  such  that  the  wavelength  component  of  550  nm  of  the  DMD  reflected  light  can  go  through 
the  double-Amici  prism  without  angular  or  spatial  deviations  from  the  optical  axis.  Other  wavelength 
components  suffer  from  different  levels  of  angular  deviations  when  propagating  through  the  double-Amici 
prism.  We  used  CDGM  glasses  H-ZF6  (refractive  index  n  =  1.7552,  Abbe  number  v d  =  24.5474)  and 
H-ZK50  (refractive  index  n  =  1.607,  Abbe  number  Vd  =  56.6572)  to  build  the  high  and  low  dispersion 
wedge  prisms.  Figure  5.11(c)  shows  the  manufactured  double-Amici  prism  and  its  mechanical  holder. 

Figure  5.12  shows  the  Zemax  spot  diagram  of  the  relay  lens/double- Amici  prism  structure  at  0  field 
angle.  In  this  spot  diagram,  we  can  see  that  due  to  the  dispersion  performance  of  the  double-Amici 
prism,  the  wavelength  components  of  450  nm  and  650  nm  are  spatially  separated  by  about  220  /im.  This 
spatial  separation  spans  over  approximately  22  CCD  pixels  and  thus  results  in  22  spectral  channels  in 
the  reconstructed  image  cubes.  The  number  of  spectral  channels  in  the  reconstructed  image  cubes  can 
be  increased  or  reduced  by  moving  the  prism  closer  to  or  further  away  from  the  relay  lens  [56].  In  our 
experimental  setup,  we  usually  align  the  optics  to  ensure  the  acquisition  of  24  spectral  images  in  the 
reconstructed  image  cubes.  From  Fig.  5.12,  we  can  also  observe  a  non-linear  dispersion  performance  of 
the  double-  Amici  prism.  We  observe  that  the  spatial  separation  between  the  wavelength  components  of 
450  nm  and  500  nm  (50  nm  difference)  is  88  /xm,  whereas  the  spatial  separation  between  the  wavelength 
components  of  600  nm  and  650  nm  (also  50  nm  difference)  is  only  about  30  fim.  We  can  see  that  this 
prism  is  more  dispersive  in  the  short  wavelength  region  of  the  visible  spectrum  than  in  the  long  wavelength 
region.  This  non-linear  dispersion  performance  was  also  discussed  in  [51]. 

We  experimentally  realized  the  DMD-SSI  in  our  lab,  which  is  shown  in  Fig.  5.13.  The  DMD  used  in 
this  setup  is  a  TI  Discovery  1100  series  DMD  product.  The  imaging  lens  is  a  Leica  COLORPLAN-P2 
projector  lens,  which  has  a  good  control  over  chromatic  aberrations.  The  focal  length  of  the  imaging 
lens  is  90  mm  and  its  f-number  is  2.5.  The  CCD  camera  is  an  AVT  Stingray  F031B  camera,  which  has 
a  pixel  dimension  of  540  x  640  and  a  pixel  pitch  of  9.9  jim.  In  this  experimental  setup,  we  use  a  2  x  2 
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Figure  5.12:  Spot  diagram  of  the  relay  lens/double- Amici  prism  structure  at  0  field  angle. 


Figure  5.13:  Experimental  setup  of  the  DMD-SSI  system. 


DMD  super-pixel  (27.36  x  27.3 6/im2)  to  represent  one  pixel  in  the  CS  measurement  pattern,  and  image 
one  DMD  super-pixel  onto  a  3  x  3  CCD  super-pixel  (29.7  x  29.7 fim2). 


5.5  Calibration  of  the  CASSI  System 

The  simulation  result  presented  in  Section  5.3  demonstrated  that  the  TwIST  algorithm  can  effectively 
reconstruct  the  original  image  cube  using  CS  measurement  results  generated  from  an  ideal  CASSI  sys¬ 
tem.  However,  in  experimental  settings,  practical  issues  such  as  noise  contaminations  and/or  optical 
aberrations  make  the  collection  of  ideal  CS  measurement  results  impossible.  In  the  following  example, 
we  present  a  simple  case  demonstrating  that  external  noise  poses  negative  impact  on  the  image  recon¬ 
struction  quality.  We  also  introduce  a  simple  calibration  process  can  be  used  to  alleviate  that  negative 
impact. 


Figure  5.14:  (a)  Partial  view  of  an  ideal  CS  measurement  pattern  (32  x  32  pixels),  (b)  A  noise- 

contaminated  CS  measurement  pattern  generated  by  adding  white  noise  to  the  pattern  shown  in  Fig. 
5.14(a). 


Figure  5.14(a)  shows  a  partial  view  (32  x  32  pixels)  of  an  ideal  CS  measurement  pattern,  which  is 
a  256  x  256  binary  random  pattern.  The  white  pixels  in  this  pattern  represent  the  binary  value  1,  and 
the  black  pixels  represent  the  binary  value  0.  In  Fig.  5.14(b),  we  added  white  noise  to  the  ideal  CS 
measurement  pattern,  which  models  the  noise  happened  on  the  DMD  mirror  plane  due  to  issues  like  dust 
scattering  and  deformations/vibrations  of  mirror  pixels  happened  in  their  tilting/resetting  process.  The 
average  intensity  of  the  white  noise  considered  in  this  simulation  is  0.53,  and  its  standard  deviation  is 
0.12.  We  used  the  ideal  and  the  noise-contaminated  patterns  to  generate  CS  measurement  results,  which 
were  then  input  into  the  TwIST  algorithm  for  reconstructing  the  original  image  cube.  For  the  sake  of 
paper  space,  we  do  not  show  the  entire  reconstructed  image  cubes.  Instead,  in  top  row  images  of  Fig. 
5.15,  we  show  RGB  images  of  the  original  (Fig.  5.15(a))  and  the  reconstructed  image  cubes  (Figs.  5.15(b) 
and  (c)).  In  bottom  row  images  of  Fig.  5.15,  we  show  optical  images  of  the  original  (Fig.  5.15(d))  and 
the  reconstructed  image  cubes  (Figs.  5.15(e)  and  (f)). 

Comparing  Fig.  5.15(b)  with  Fig.  5.15  (c)  and  comparing  Fig.  5.15(e)  with  Fig.  5.15(f),  we  can 
see  that  the  quality  of  the  reconstructed  RGB  and  optical  images  is  apparently  affected  by  the  noise 
considered  in  the  CS  measurement  pattern,  in  the  sense  that  Figs.  5.15(c)  and  (f)  look  more  blurry  and 
noisy  than  Figs.  5.15(b)  and  (e).  The  PSNR  value  of  the  optical  image  reconstructed  using  the  noisy 
CS  measurement  pattern  is  15.69  dB,  which  is  9.6%  lower  than  the  PSNR  value  of  the  optical  image 
reconstructed  using  the  ideal  CS  measurement  pattern  (17.36  dB).  To  mitigate  the  quality  degradation 
caused  by  noise  contaminations,  when  solving  the  image  reconstruction  problem,  it  is  helpful  to  consider 
the  same  amount  of  noise  that  happened  in  the  CS  measurement  process  in  the  system  transfer  matrix  H. 

In  Fig.  5.16(a),  we  show  the  RGB  image  reconstructed  when  the  same  amount  of  noise  that  happened 
in  the  CS  measurement  process  was  considered  in  the  system  transfer  matrix.  Compared  to  Fig.  5.16(b), 
which  is  the  RGB  image  reconstructed  without  considering  the  noise  in  the  system  transfer  matrix  ,  we  can 
see  that  Fig.  5.16(a)  looks  less  noisy  and  the  color  contrast  between  the  red  and  green  tomatoes  is  more 
visually  apparent.  Figure  5.16(c)  shows  the  optical  image  reconstructed  when  the  noise  is  considered  in 
the  system  transfer  matrix  and  Fig.  5.16(d)  shows  the  optical  image  reconstructed  with  the  ideal  system 
transfer  matrix.  We  calculated  the  PSNR  values  of  those  reconstructed  optical  images  and  noted  that,  by 
considering  the  noise  in  the  system  transfer  matrix,  the  PSNR  value  of  the  reconstructed  optical  image 
was  improved  from  15.69  dB  (Fig.  5.16(d))  to  17.07  dB  (Fig.  5.16(c)). 

In  the  above  simulation,  we  only  considered  one  kind  of  noise  contaminations  in  the  DMD-SSI  system. 
In  experimental  settings,  other  kinds  of  noise  contaminations,  such  as  mechanical  vibrations,  stray  light 
in  the  lab  environment,  and  electronic  noise  of  the  camera,  cause  additional  errors  to  the  CS  measurement 
results.  To  fully  account  for  those  non-ideal  factors,  a  calibration  process  was  proposed  in  [51,53],  which 
involves  the  acquisition  of  a  series  of  monochromic  images  for  the  CS  measurement  pattern.  Each  of 
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Figure  5.15:  (a)  RGB  image  of  the  original  image  cube,  (b)  RGB  image  of  the  reconstructed  image  cube 
when  the  ideal  random  pattern  was  used  to  implement  the  CS  measurement  process,  (c)  RGB  image 
of  the  reconstructed  image  cube  when  the  noise-contaminated  CS  measurement  pattern  was  used,  (d) 
Optical  image  of  the  original  image  cube,  (e)  Optical  image  of  the  reconstructed  image  cube  when  the 
ideal  CS  measurement  pattern  was  used,  (f)  Optical  image  of  the  reconstructed  image  cube  when  the 
noise-contaminated  CS  measurement  pattern  was  used. 


the  monochromic  images  contains  not  only  the  spatial  information  of  the  implemented  CS  measurement 
pattern,  but  also  all  the  wavelength-dependent  noises  and  aberrations.  The  acquired  monochromic  images 
are  then  used  to  build  a  practical  system  transfer  matrix  for  the  experimental  setup,  such  that  all  the 
noises  and  aberrations  present  in  the  experiment  setup  can  be  considered  in  the  image  reconstruction 
process.  Figure  5.17  shows  the  calibration  setup  for  our  DMD-SSI  system,  which  is  composed  of  a  wide¬ 
band  Xenon  lamp  (Newport,  stock  number:  66477)  and  a  visible  monochromator  (Oriel,  model  number: 
77200). 

In  the  calibration  process,  monochromic  emissions  generated  from  the  monochromator  are  used  to 
illuminate  the  DMD  mirror  plane  and  corresponding  monochromic  images  of  the  implemented  DMD 
pattern  can  be  captured  by  the  CCD  camera  installed  in  the  DMD-SSI  system.  Figure  5.18(a)  and  (b) 
show  an  ideal  CS  measurement  pattern  and  one  of  its  monochromic  images  captured  at  the  wavelength 
of  612  nm.  In  this  figure,  we  can  see  that  the  monochromic  image  is  contaminated  with  random  noises 
and  optical  blurs.  Also,  the  intensity  contrast  between  the  on  and  off  pixels  is  lower  in  the  monochromic 
image  than  in  the  ideal  CS  measurement  pattern.  All  of  those  non-ideal  factors  will  be  considered  in  the 
image  reconstruction  process  by  the  practical  system  transfer  matrix. 

In  the  calibration  process,  we  capture  monochromic  images  every  1  nm  from  450  nm  to  685  nm. 
Therefore,  a  total  number  of  236  monochromic  images  are  captured.  Those  monochromic  images  are 
spatially  separated,  due  to  the  dispersion  performance  of  the  prism.  Not  all  of  those  236  monochromic 
images  are  used  to  build  the  transfer  matrix  for  the  experimental  setup.  This  is  because  the  CCD  camera 
cannot  accurately  resolve  spatial  separations  smaller  than  1  CCD  pixel.  Therefore,  monochromic  images 
that  have  spatial  separations  smaller  than  1  CCD  pixel  are  considered  as  identical  in  this  calibration 
process.  From  those  236  monochromic  images  captured  in  the  calibration  process,  we  found  24  non- 
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(a)  (b) 


Figure  5.16:  RGB  and  optical  images  reconstructed  using  the  noise-contaminated  CS  measurement  pat¬ 
tern.  In  (a)  and  (c),  the  same  amount  of  noise  happened  in  the  CS  measurement  process  was  considered 
in  the  system  transfer  matrix.  In  (b)  and  (d),  the  noise  was  not  considered  in  the  system  transfer  matrix. 


Figure  5.17:  Calibration  setup  for  the  DMD-SSI  system. 


identical  images  and  used  them  to  build  the  practical  system  transfer  matrix.  We  say  those  24  non¬ 
identical  monochromic  images  form  a  calibration  cube  for  the  experimental  setup.  Figure  5.19(a)  shows 
the  wavelength  information  of  the  24  non-identical  monochromic  images.  Figure  5.19(b)  shows  the  first 
and  the  last  monochromic  images  included  in  the  calibration  cube,  which  were  captured  at  the  wavelengths 
of  453  nm  and  671  nm,  respectively.  Fake  colors  are  added  to  those  two  monochromic  images  to  enhance 
their  visual  perception.  We  can  see  that  those  two  monochromic  images  are  spatially  separated  by  23 
CCD  pixels  in  the  horizontal  direction  of  the  CCD  camera.  Please  note  that  the  calibration  process  is 
only  useful  to  one  experimental  setup.  If  the  optical/mechanical  architecture  of  the  DMD-SSI  system  is 
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Figure  5.18:  (a)  An  ideal  CS  measurement  pattern  (128  x  128).  (b)  A  monochromic  image  of  the  CS 
measurement  pattern  captured  by  the  DMD-SSI  system  at  612  nm. 


modified,  such  as  a  new  lens  is  installed  or  the  position  of  the  prism  is  moved,  a  new  calibration  process 
needs  to  be  implemented  for  the  modified  setup. 
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Figure  5.19:  (a)  Wavelength  information  of  non-identical  monochromic  images  captured  in  the  calibration 
process,  (b)  Monochromic  images  captured  at  453  nm  (red  pattern)  and  671  nm  (blue  pattern). 


5.6  Experimental  Results 

Figure  5.20  shows  an  image  reconstruction  result  obtained  from  the  experimental  DMD-SSI  setup.  In 
this  case,  the  binary  random  pattern  shown  in  Fig.  5.18  was  used  as  the  CS  measurement  pattern.  Figure 
5.20(a)  shows  a  CCD  image  of  the  imaging  target,  which  is  composed  of  a  green  and  a  red  toy  car.  Figure 
5.20(b)  shows  the  CS  measurement  result  generated  from  the  DMD-SSI  setup.  Figure  5.20(c)  shows  the 
reconstructed  image  cube  of  the  imaging  target.  We  can  see  that  the  green  car  was  reconstructed  with 
strong  intensities  in  green  spectral  channels  (529  nm  to  580  nm),  whereas  the  red  car  was  reconstructed 
with  strong  intensities  in  red  spectral  channels  (605  nm  to  636  nm). 

In  Fig.  5.21,  we  show  another  set  of  image  reconstruction  result  which  was  generated  with  a  different 
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Figure  5.20:  (a)  CCD  image  of  the  imaging  target,  (b)  CS  measurement  result  acquired  by  the  DMD-SSI 
system  for  the  imaging  target,  (c)  Reconstructed  spatial/spectral  image  cube  (fake  colors  were  added  to 
enhance  the  visual  perception  of  those  spectral  images). 


imaging  target.  Also,  a  different  CS  measurement  pattern  was  used  to  realize  the  CS  measurement 
process.  A  new  calibration  process  was  implemented  prior  to  this  image  reconstruction  experiment. 
Figure  5.21(a)  shows  a  monochromic  image  of  the  CS  measurement  pattern  used  in  this  experiment, 
which  was  captured  at  the  wavelength  of  587  nm.  Figure  5.21(b)  shows  a  CCD  image  of  the  imaging 
target  used  in  this  experiment,  which  is  composed  of  three  colorful  plastic  stripes.  Figure  5.21(c)  shows 
the  CS  measurement  result  captured  from  the  experimental  DMD-SSI  setup.  Figure  5.21(d)  shows  the 
reconstructed  spatial/spectral  image  cube  of  the  imaging  target.  In  this  result,  we  can  see  that  the  blue 
color  stripe  was  reconstructed  with  strong  intensities  in  spectral  channels  from  460  nm  to  498  nm,  whereas 
the  green  color  stripe  was  reconstructed  with  strong  intensities  in  spectral  channels  from  517  nm  to  576 
nm  and  the  red  color  stripe  was  reconstructed  with  strong  intensities  in  spectral  channels  from  601  nm 
to  663  nm.  This  reconstruction  result  is  in  accordance  with  our  visual  perception  to  the  blue,  green,  and 
red  color  stripes  in  the  imaging  target. 

5.7  Multi-shot  CS  Measurements 

Multi-Shot  CASSI  (MS-CASSI)  systems  were  developed  to  enhance  the  quality  of  image  cubes  recon¬ 
structed  by  conventional  Single-Shot  CASSI  (SS-CASSI)  systems.  In  MS-CASSI  systems,  the  photomask 
is  mounted  on  a  piezo  stage.  By  horizontally  and/or  vertically  shifting  the  photomask  using  the  stage, 
different  regions  of  the  random  pattern  coated  on  the  photomask  can  be  used  to  impose  intensity  modula¬ 
tions  on  the  original  image  cube.  MS-CS  measurement  processes  can  be  realized  by  using  multiple  shifts 
of  the  piezo-stage.  In  [52],  the  authors  demonstrated  the  quality  enhancement  brought  by  the  MS-CASSI 
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Figure  5.21:  (a)  Monochromic  image  of  the  CS  measurement  pattern  captured  at  587  nm.  (b)  CCD 
image  of  the  imaging  target,  (c)  CS  measurement  result,  (d)  Reconstructed  spatial/spectral  image  cube 
of  the  imaging  target  (fake  colors  were  added  to  enhance  the  visual  perception). 


system.  DMD-SSI  systems  can  also  be  used  to  realize  MS-CS  measurement  processes  [55].  Compared  to 
piezo  stage-based  MS-CASSI  systems,  DMD-SSI  systems  realize  MS-CS  measurement  processes  without 
mechanical  motions,  and  thus,  they  are  free  from  registration  errors  happened  in  the  motion  process. 
More  importantly,  DMD-SSI  systems  has  a  much  higher  flexibility  in  implementing  CS  measurement 
patterns,  in  the  sense  that  in  addition  to  the  shifted  random  patterns  that  the  piezo  stage-based  system 
can  only  provide,  DMD-SSI  systems  can  use  patterns  of  independent  designs  to  realize  MS-CS  measure¬ 
ment  processes. 

We  implemented  MS-CS  measurement  processes  with  our  DMD-SSI  setup.  Figure  5.22(a)  shows  the 
imaging  target  used  in  this  experiment,  which  is  a  red  chili  pepper  with  a  green  stem.  Figure  5.22(b)  shows 
one  of  the  CS  measurement  patterns  used  in  this  experiment.  Figure  5.22(c)  shows  the  CS  measurement 
result  generated  when  the  CS  measurement  pattern  shown  in  Fig.  5.22(b)  was  implemented  with  the 
DMD.  Top  row  images  in  Fig.  5.23  show  spectral  images  (channels  13-24  of  the  reconstructed  image 
cube)  reconstructed  with  a  Single-Shot  CS  (SS-CS)  measurement  process.  Center  row  images  in  Fig. 
5.23  show  spectral  images  reconstructed  with  a  6-shot  MS-CS  measurement  process  implemented  with 
shifted  random  patterns  and  bottom  row  images  show  spectral  images  reconstructed  with  a  6-shot  MS-CS 
measurement  process  implemented  with  independent  random  patterns. 

In  this  figure,  we  can  see  that  both  the  SS-CS  and  the  MS-CS  measurement  processes  can  result  in 
meaningful  image  reconstruction  results.  We  can  see  that  the  green  stem  portion  of  the  pepper  target 
was  reconstructed  with  strong  intensities  in  green  spectral  channels,  which  are  from  538  nm  to  581  nm 
(channels  14-18),  and  the  red  body  portion  of  the  pepper  target  was  reconstructed  with  strong  intensities 
in  red  spectral  channels,  which  are  from  606  nm  to  670  nm  (channels  20-24).  However,  we  noted  that 
the  spectral  images  reconstructed  with  MS-CS  measurement  processes  have  much  better  quality  than  the 
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Figure  5.22:  (a)  CCD  image  of  the  imaging  target,  (b)  One  of  the  CS  measurement  patterns  used  in  the 
MS-CS  measurement  experiment,  (c)  CS  measurement  result  generated  from  our  experimental  setup, 
when  the  pattern  shown  in  Fig.  5.22(b)  was  implemented  with  the  DMD. 


optical  image  (529  nm)  (538  nin)  (548  nm)  (558  nm)  (569  nm)  (581  nm)  (593  nm)  (606  nm)  (621  nm)  (637  nm)  (653nm)  (670nm) 

Figure  5.23:  Top  row  images:  spectral  images  (channels  13-24  in  the  reconstructed  image  cube)  recon¬ 
structed  with  a  SS-CS  measurement  process.  Center  row  images:  spectral  images  reconstructed  with  a 
6-shot  MS-CS  measurement  process,  implemented  with  shifted  random  patterns.  Bottom  row  images: 
spectral  images  reconstructed  with  a  6-shot  MS-CS  measurement  process,  implemented  with  independent 
random  patterns. 


spectral  images  reconstructed  with  the  SS-CS  measurement  process,  in  the  sense  that  they  look  less  noisy 
and  blurry.  Also,  in  spectral  images  reconstructed  with  MS-CS  measurement  processes,  the  intensity 
contrast  between  the  body /stem  portions  of  the  pepper  target  and  the  dark  background  is  more  visually 
apparent. 

In  Fig.  5.24,  we  show  three  colorful  images  generated  by  using  channel  16  and  channel  21  in  the 
reconstructed  image  cubes  as  the  green  and  red  color  components  in  the  RGB  data  format.  We  can  see 
that  those  colorful  images  look  similar  to  the  CCD  image  of  the  imaging  target  shown  in  Fig.  5.22(a). 
Also  in  this  figure,  we  can  see  that  colorful  images  reconstructed  with  MS-CS  measurement  processes 
look  less  noisy  than  the  one  reconstructed  with  the  SS-CS  measurement  process.  They  also  have  a  higher 
color  contrast  between  the  stem  and  the  body  portions  of  the  pepper  target.  We  also  evaluated  the 
intensity  variations  of  the  reconstructed  image  cubes  at  point- 1  and  point-2  on  the  imaging  scene,  which 
are  located  on  the  body  and  stem  portions  of  the  pepper  target,  respectively.  The  locations  of  point- 1  and 
point-2  are  demonstrated  in  Fig.  5.22(a).  Those  intensity  variations  can  be  considered  as  reconstructed 
spectral  curves  on  those  two  pixel  locations.  We  also  used  a  commercial  spectrometer  (Avantes,  model 
number:  AvaSpec-1024)  to  generate  reference  spectral  curves  at  those  two  locations.  Figure  5.25  shows 
the  reconstructed  (dotted  lines)  and  the  reference  spectral  curves  (solid  lines). 
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(a)  (b)  (c) 


Figure  5.24:  Colorful  images  synthesized  using  spectral  channels  16  and  21  in  the  reconstructed  image 
cubes  as  the  green  and  red  color  components  in  the  RGB  data  format,  (a)  Colorful  image  generated 
with  the  SS-CS  measurement  process,  (b)  Colorful  image  generated  with  a  6-shot  MS-CS  measurement 
process,  implemented  with  shifted  random  patterns,  (c)  Colorful  image  generated  with  a  6-shot  MS-CS 
measurement  process,  implemented  with  independent  random  patterns. 
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Figure  5.25:  (a)  Reconstructed  and  reference  spectral  curves  measured  at  point- 1  on  the  pepper  target, 
(b)  Reconstructed  and  reference  spectral  curves  measured  at  point-2  on  the  pepper  target. 


In  Fig.  5.25,  we  can  see  that  the  reference  spectral  curves  (solid  curves)  have  an  intensity  peak  at 
the  spectral  channel  of  621  nm  for  point- 1  and  an  intensity  peak  at  the  spectral  channel  of  548  nm  for 
point-2.  We  call  those  two  spectral  channels  the  fingerprint  spectral  channels  of  the  pepper  target  in  the 
body  and  stem  portions.  In  spectral  curves  reconstructed  with  the  SS-CS  measurement  process  (dotted 
line  in  dark  yellow  color),  we  can  see  the  intensity  peaks  miss  the  ones  in  the  reference  spectral  curves, 
which  implies  the  reconstruction  quality  is  not  ideal.  In  this  case,  the  MSE  values  of  the  deviations 
from  the  reconstructed  spectral  curves  to  the  reference  curves  are  3.3%  at  point-1  and  2.9%  at  point-2, 
respectively.  For  spectral  curves  reconstructed  with  MS-CS  measurement  processes  (dotted  curves  in  blue 
and  pink  colors),  we  can  see  that  their  intensity  peaks  match  the  ones  in  the  reference  spectral  curves. 
Also,  MSEs  of  the  deviations  from  the  reconstructed  spectral  curves  to  the  reference  curves  dropped  to 
1.2%  and  1.4%  for  point-1,  and  to  2.4%  and  0.4%  for  point-2,  when  shifted  and  independent  random 
patterns  were  used  to  implement  the  MS-CS  measurement  processes,  respectively. 
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Chapter  6 


Communications  System 


The  objective  of  this  task  is  to  develop  a  robust,  low-delay,  low-complexity,  communication  system  for 
the  transmission  of  multispectral  measurements.  In  order  to  achieve  robustness  against  changes  in  the 
channel  quality,  and  to  reduce  the  system  complexity  with  respect  to  standard  communications  systems, 
we  have  developed  a  discrete-time,  analog-processing  system  that  completely  skips  the  digital  domain  and 
achieves  a  performance  close  to  the  theoretical  limits.  The  key  idea  at  the  encoder  site  is  to  perform  simple 
non-linear  transformations  on  the  available  measurements.  These  transformations,  denoted  as  M()  and 
T()  in  Fig.  1,  depend  on  the  source/channel  statistics  and  on  the  desired  transmission  rate.  The  resulting 
values  are  directly  transmitted  though  the  communications  channel  using  very  simple  PAM  signaling.  ML 
decoding  at  the  receiver  is  straightforward,  and  leads  to  excellent  performance.  The  research  developed 
in  this  project  has  focused  on  the  development  of  the  appropriate  non-linear  mappings,  and  on  the 
implementation  of  testbeds  to  study  the  performance  of  the  proposed  system  in  practical  scenarios. 


Encoder 


a  ix^) 


Figure  6.1:  Discrete-time  continuous-amplitude  communications  system 


6.1  Development  of  non-linear  mappings 

A  key  piece  of  our  work  is  the  development  of  low-complex  non-linear  mappings  that  match  the  statistics 
of  the  measurements  to  the  channel,  leading  to  performance  close  to  the  theoretical  limits.  We  denote  the 
number  of  measurements  to  be  transmitted  as  N,  and  K  denotes  the  total  number  of  allowed  channel  uses. 
Then,  the  theoretical  limit  establishing  the  minimum  amount  of  distortion  (or  equivalently,  the  maximum 
Signal  to  Distortion  Ratio,  SDR)  that  can  be  achieved  for  a  certain  channel  quality  is  a  function  of  the 
bandwidth  reduction  factor  N/K  and  the  Rate  Distortion  Function,  R(D).  This  theoretical  limit  holds 
for  any  communication  system,  and  our  goal  is  to  get  as  close  as  possible  to  the  limit  with  the  lowest 
possible  complexity. 

We  started  our  research  with  the  development  of  non-linear  mappings  for  the  case  of  independent  and 
identically  distributed  measurements  transmitted  through  RF  wireless  channels.  We  also  implemented 
quasi-optimal  decoding  algorithms  with  very  reduced  complexity,  close  to  that  of  Maximum  Likelihood, 
ML.  Next,  we  investigated  mappings  that  exploit  the  correlation  existing  between  the  measurements 
to  achieve  an  improved  performance.  Finally,  we  also  developed  mappings  for  the  case  in  which  the 


measurements  are  transmitted  optically,  which  presents  new  research  problems  as  the  characteristics  of 
optical  wireless  channels,  and  thus  the  mapping  designs,  are  very  different  from  that  of  RF  channels. 
Through  the  process,  we  also  discovered  a  technique  to  simplify  even  more  the  encoding  process  when 
the  bandwidth  reduction  factor  is  2,  decreasing  the  computational  complexity  of  the  system  by  an  order 
of  magnitude.  Interestingly,  the  resulting  performance  in  this  case  is  very  similar  to  that  obtained  with 
the  original  mappings. 

6.2  Testbeds  implementation 

In  addition  to  the  theoretical  work,  we  have  implemented  two  testbeds  to  study  the  performance  of  the 
proposed  system  in  practical  scenarios.  The  first  testbed,  which  will  be  described  in  more  detail  in  the 
sequel,  uses  Software  Defined  Radios  (SDRs)  for  the  electromagnetic  transmission  of  the  measurements 
through  wireless  channels,  while  the  second  testbed  performs  wireless  optical  communications.  The  RF 
testbed  is  composed  of  two  Universal  Software  Radio  Peripherals  (USRP2)  (see  Fig.  2).  The  transmit  and 
receive  antennas  are  separated  2  meters,  and  thus  the  wireless  channel  is  a  quasi-static  fading  channel, 
which  for  the  transmission  of  each  frame  behaves  as  an  Additive  White  Gaussian  Noise  (AWGN)  channel. 
The  system  has  a  frequency  of  operation  in  the  ISM  band,  fc=2.41  GHz,  and  the  baseband  sampling 
frequency  is  fs=195  samples/s.  We  apply  a  squared  root  raised  cosine  pulse-shaping  filter  with  a  roll-off 
factor  of  20%  and  the  symbol  period  is  Ts=20  samples/s.  As  explained  before,  we  used  the  very  simple 
PAM  modulation  and  transmit  data  just  using  the  in-phase  component. 


(a)  Transmitter  (b)  Receiver 


Figure  6.2:  RF  testbed 

In  order  to  perform  frequency  synchronization,  the  testbed  used  a  global  GPS  Disciplined  Oscillator 
shared  between  the  transmitter  and  the  receiver.  The  main  goal  of  this  procedure  is  to  align  the  frequency 
of  the  carrier  at  the  receiver  with  the  frequency  of  the  carrier  at  the  transmitter.  The  result  is  a 
coarse-grained  frequency  synchronization  that  avoids  the  continue  rotation  of  the  received  symbols  at 
the  receiver.  Although  this  technique  solved  in  general  the  problem  of  frequency  mismatch,  there  was  a 
remaining  impairment  related  to  a  constant  carrier  phase  error.  We  addressed  this  impairment  using  SISO 
frequency  estimator  algorithms  (Kays  estimator)  that  leads  to  good  results  in  the  testbed  implementation. 
In  terms  of  frame  synchronization,  we  divide  the  data  in  frames  of  fixed  size  (10,000  symbols),  and  at  the 
beginning  of  each  frame  we  add  a  known  pilot  sequence  of  101  symbols  generated  by  a  BPSK  modulation. 
At  the  beginning  of  each  received  frame,  we  perform  a  cross  correlation  procedure  between  the  received 
frame  and  the  expected  pilot  sequence  in  order  to  find  the  start  of  the  sequence. 


6.3  Results 

Our  results  focus  on  the  transmission  of  the  multispectral  measurements  through  the  RF  testbed.  As 
explained  before,  the  development  of  appropriate  mappings  depends  on  the  statistics  of  the  measurements, 
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(a)  Original  measurements 


xio4 


(b)  Their  distribution 


Figure  6.3:  Notice  the  good  match  between  the  measurements  and  the  Gaussian  distribution 


which,  as  shown  in  Fig.  3,  can  be  accurately  approximated  by  a  Gaussian  distribution.  We  consider  two 
bandwidth  reduction  factors:  N/K=l,  or  direct  transmission,  and  N/K=2. 

In  the  case  of  N/K=l,  it  is  well  known  that  direct  transmission  through  an  AWGN  channel  is  optimal 
provided  that  the  measurements  are  iid  and  Gaussian.  Therefore,  since  the  wireless  channel  can  be 
approximated  as  AWGN,  the  performance  of  the  proposed  system  should  approach  the  theoretical  limit 
for  Gaussian  sources,  which  can  be  easily  calculated  using  information  theoretical  arguments.  Indeed, 
Fig.  4 (left)  shows  that  simulation  results  with  MMSE  decoding  perfectly  fit  the  theoretical  limit,  while, 
as  expected,  ML  decoding  results  in  performance  lost  for  low  signal  to  noise  ratios  in  the  channel  (CSNR). 
Notice  the  excellent  match  between  simulation  results  and  the  testbed  evaluation.  Similar  results  can  be 
seen  in  Fig.  4(right)  for  the  case  in  which  the  bandwidth  reduction  factor  is  N/K=2. 


Figure  6.4:  Theoretical  limit  and  performance  of  the  proposed  system  for  the  transmissiosn  of  the  original 
measurements  through  an  RF  wireless  channel.  The  system  is  evaluated  through  simulations  and  through 
the  testbed. 
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